古典迭代法中的Jacobi迭代法和SOR迭代法求解Hilbert矩阵方程组.pdf
优质文本燕山大学课课 程程 设设 计计 说说 明明 书书题目：直接法求解题目：直接法求解 HilbertHilbert 矩阵方程组矩阵方程组学院系学院系：理学院理学院年级专业：年级专业：学学号：号：学生姓名：学生姓名：指导教师：指导教师：教师职称：教师职称：教授教授1/37优质文本燕山大学课程设计论文任务书燕山大学课程设计论文任务书院系：理学院基层教学单位：学生姓学 号设计题古典迭代法求解 Hilbert 矩阵方程组目设计技术参数专业 班级名分别利用古典迭代法中的 Jacobi 迭代法和 SOR 迭代法，求解 Hilbert 矩阵方程组，并及准确解进行比拟，得出结论，然后针对不同的 n 值重复进行计算并比拟所得的结果。写出相应的 Matlab 程序，执行程序得出的结果。通过结果比照分析并进行理论分析。设计要求四周的时间掌握 Hilbert 矩阵的性质、Hilbert 矩阵方程组工作量的特点及其相关的求解方法，研究古典迭代法中的Jacobi迭代法和SOR迭代法并用这两种方法求解Hilbert矩阵方程组，分析这两种方法对解 Hilbert 矩阵方程组的效果如何；熟悉使用 Matlab 软件，上机求解给定题目。2/37优质文本第一周：查阅资料，根本了解 Jacobi 迭代法和 SOR 迭代法的根本原理、解题步骤等，主要对 Jacobi 迭代法进行分析；工第二周至第三周：根据其他组员编写的 Jacobi 迭代法的作Matlab 程序进行具体题目的计算并及其他同学负责的局部计划进行综合分析；第四周：整理和撰写报告，进行辩论。1谢进，李大美.MATLAB 及计算方法实验M.武汉大学出版社，2016.2周品，何正风.MATLAB 数值分析M.机械工业出版社，参考资料2016.13林雪松，周婧，林德新.MATLAB7.0 应用集锦M.机械工业出版社，20064(美)John H.Mathews,Kurtis D.Frink 数值方法M.电子工业出版社，20055李庆扬，王能超，易大义.数值分析M.华中科技大学出版社，2006.1指导教师基层教学单位主签字任签字说明：此表一式四份，学生、指导教师、基层教学单位、系部各一份。3/37优质文本2012 年 01 月 09 日燕山大学课程设计评审意见表指导教师评语：成绩：指导教师：年月日辩论小组评语：成绩：评阅人：4/37优质文本年月日课程设计总成绩：辩论小组成员签字：年月日5/37优质文本目目录录摘要2一、一、课题背景简介3二、数学模型描述3三、算法原理4 (一)Hilbert 矩阵的来源及其正定性.4 (二)迭代法的根本思想.51Jacobi 迭代法.62SOR 方法 7四、算法程序8(一)Jacobi 迭代法8(二)SOR 方法 9五、求解过程10一第1题10二第2题12三第3题18六、心得体会20七、参考文献201/37优质文本摘要摘要论文旨在阐述求解 Hilbert 矩阵方程方法中的古典迭代法应用于求解特殊矩阵的过程,其中主要包括Jacobi迭代法和SOR迭代法。文中从方法原理、算法的根本应用到理论的归纳及扩展。通过这一例题，可进一步领会数值分析的实际应用。同时为了考察读者对古典迭代法的学习灵活度，例题的解答需要用到数值方法在MATLAB 中的实现，包括数值方法在 MATLAB 中的函数实现，MATLAB编程和 MATLAB 绘图，使读者在上机练习中加深对数值分析算法原理的理解，通过对算法的思考和理论分析，掌握MATLAB 的使用。关键词关键词：Hilbert 矩阵古典迭代法 Jacobi 迭代法 SOR 迭代法2/37优质文本一、课题背景简介随着电子计算机的出现和迅速开展，在各门自然科学和工程技术科学的开展中，科学计算已经成为平行于理论分析和科学实验的第三种科学手段。数值计算是科学计算中的一个不可缺少的环节，而在数值计算中，一类很重要的问题就是线性方程组的求解。另外，数学和物理以及力学等学科和工程技术中许多问题的最终解决都归结为解一个或一些大型系数矩阵的线性方程组。所以，大型线性方程组的求解是大规模科学及计算的核心，许多作者都对此作了研究。对线性方程组Axb的求解，主要有直接求解法和迭代法求解。对于阶数不太高的线性方程组，用直接法比拟方便，高斯消元法和克兰姆法那么是直接解法里面最重要的解法。但是，随着科学技术的飞速开展，需要求解的问题的规模越来越大，迭代法已取代直接法成为求解大型线性组的最重要的一类方法。此时，迭代格式3/37优质文本的收敛性和收敛速度成为一个很重要的问题，成为人们关注的焦点。不收敛的格式当然不能用，虽然收敛但是收敛的很慢的格式，不仅是人工和机器的时间比拟浪费，而且还不一定能解出结果，实际应用价值太小。通常用的迭代法有 Jacobi，Gauss-Seidel 等古典迭代法，还有SOR，AOR 以及 SAOR，SSOR 等迭代法。这些迭代法的提出对大型线性方程组的求解提供了一条快速有效的途径。但是这些迭代法相对来说使用条件很苛刻，或者很繁琐。此外当迭代矩阵的谱半径比拟大，尤其接近 1 时，迭代速度将会很慢，极大的影响了方程组求解的时效性，从而使它们的应用受到了一定的限制。MATLAB 是矩阵实验室Matrix Laboratory的简称，是美国MathWorks 公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括 MATLAB 和 Simulink 两大局部。MATLAB 和 Mathematica、Maple 并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB 可以进行矩阵运算、绘制函数和数据、实现算法、创立用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理及通讯、图像处理、信号检测、金融建模设计及分析等领域。MATLAB 的根本数据单位是矩阵，它的指令表达式及数学、工程中常用的形式十分相似，故用 MATLAB 来解算问题要比用 C，4/37优质文本FORTRAN 等语言完成相同的事情简捷得多，并且 mathwork 也吸收了像 Maple 等软件的优点,使 MATLAB 成为一个强大的数学软件。在新的版本中也参加了对 C，FORTRAN，C+，JAVA 的支持。可以直接调用,用户也可以将自己编写的实用程序导入到 MATLAB 函数库中方便自己以后调用，此外许多的 MATLAB 爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。MATLAB 的优势：1友好的工作平台和编程环境；2简单易用的程序语言；3强大的科学计算机数据处理能力；4出色的图形处理功能；5应用广泛的模块集合工具箱；6 实用的程序接口和发布平台；7应用软件开发包括用户界面二、数学模型描述：二、数学模型描述：古典的迭代法求解 Hilbert 矩阵方程组考虑方程组Hx b,其中系数矩阵为：H (hij)Rnn,hij1,i,j 1,2,ni j 1假设由准确解x(1,1,1)T Rn确定向量b。1选择 n=6，分别用 Jacobi 迭代法和 SOR 迭代法1,1.25,1.5等求解。将计算结果及准确解进行比拟。2逐步增大n(n=8,10,),重复1的计算。比拟结果，结果5/37优质文本说明什么？3对于计算过程中发现的问题，尝试用其他方法加以解决。三、算法原理三、算法原理:(一)Hilbert 矩阵的来源及其正定性设f(x)C0,1,求 n 次多项式P(x)a0 a1x a2x2 anxn使得L P(x)f(x)2dx01取最小值。其中，P(x)称为函数f(x)的最正确平方逼近多项式。由极值必要条件，对L(a0,a1,an)ajxj f(x)2dx0j01n中的变量求导数，得nn111Liji j 2x ajx f(x)dx 2ajxdx 2xif(x)dx000aij0j0令其为零，得方程组11a xif(x)dx(i 0,1,n)j0j0i j 1n令bi0 xif(x)dx，将方程组写为矩阵形式1/211/21/31/(n1)1/(n 2)1/(n1)a0b0ab1/(n 2)11 1/(2n1)anbn6/371优质文本方程组的系数矩阵就是著名的 Hilbert 矩阵。多项式可表示为内积形式1 xP(x)a0a1 an nx 所以1 a0 xaa1 an1 x xn1 nxanP(x)2a0积分，得1/211/21/3a1 an1/(n 1)1/(n 2)1/(n 1)a0a1/(n 2)1 1/(2n1)an2P(x)dx a001由此可知，Hilbert 矩阵是对称正定矩阵。另外，Hilbert 矩阵是著名的病态矩阵。(二)迭代法的根本思想对于大型方程组Ax b1其中ARnn，x，bRn。假设 A 非奇异，方程组 有唯一解Tx(x1,x2,xn)。将Ax b变形为等价的方程组7/37优质文本x Bx f2由此建立迭代公式x(k1)Bx(k)f,k 0,1,2,3给定初始向量x(0)，按此公式计算得近似解向量序列x(k)，称此方法为迭代法。假设对任意x(0)，当迭代次数无限增加时，序列x(k)都x(k)x，也就是说limx(k)x，i 1,2,n。有相同的极限x，即limkk显然有x Bx f4这是说迭代公式是收敛的，否那么说发散的，称迭代格式3中的矩阵 B 为迭代矩阵，对于不同的迭代矩阵得到不同的迭代格式。1Jacobi 迭代法设有方程组 (i 1,2,n),记作Ax b,1A 为非奇异矩阵aij 0(i 1,2,n)。将 A 分解为A D LU,其中a11D ann 0a21L a31an10a320an2an,n10a22，8/37优质文本0a120U a1na23a2n0an1,n0a13将式1的第i(i 1,2,n)个方程用aij去除再移项，得到等价方程n1组xi(biaijxj)(i 1,2,n)，2aiij1ji简记作x B0 x f,其中B0 I D1A D1(L U),f D1b。对方程组2应用迭代法，得到解式1的 Jacobi 迭代公式(0)(o)Tx(0)(x10,x2,xn)(初始向量)，n1(k1)3xi(biaijx(jk)，aiii1ji(k)(k)T,xn)为第 k 次迭代向量。设x已经算出，由其中x(k)(x1k,x2(k)式(3)可计算下一次迭代向量x(k1)k 0,1,2,;i 1,2,n。显然迭代公式3的矩阵形式为其中B0称为 Jacobi 迭代矩阵。2SOR 方法设有方程组Ax b,1其中AR为非奇异矩阵，且设aii 0i 1,2,n，分解 A 为A D LU。29/37nn优质文本设第 k 次迭代向量x，及第 k+1 次迭代向量xj 1,2,i 1,要求计算分量xi(k1)。首先用 Gauss-Seidel 迭代法定义辅助量i11(k1)xi(biaijx(jk1)aiij1(k)(k1)的分量x(jk1)ji1(k)a xijj)(i 1,2,n)，3nxi(k1)某个平均值及加权平均再把xi(k1)取为xi(k)及，即xi(k!)(1)xi(k)xi(k!)xi(k)(xi(k!)xi(k)4用式(3)代入式4即得到解方程组Ax b的逐次超松弛迭代公式i1n(k!)(k)(k!)xi(biaijxjaijx(jk),xi5aiij1jix(k)(x(k),x(k),x(k)T(k 0,1,;i 1,2,n),12n其中称为松弛因子，或写为xi(k!)xi(k)xi(k 0,1,;i 1,2,n),i1nx(b a x(k1)a x(k).6iiijjijjaj1jiii显然，当1时，解式1的 SOR 方法就是 Gauss-Seidel迭代法。当1时，称式5为低松弛法，当1时称式5为超松弛法。SOR 迭代公式的矩阵形式，迭代公式5也可写为aiix(k1)i(1)aiix(k)i(biaijxj1i1(k1)jji1an(k)xijj)(i 1,2,n)7用式A D LU，那么Dx(k1)(b Lx(k1)Ux(k)(1)Dx(k)，10/37优质文本即(D L)x(k1)(1)D U)x(k)b.显然对于任何一个值，(D L)非奇异由设aii 0;i 1,2,n，于是x(k!)(D L)1(1)D Ux(k)(D L)1b.也就是说，解式1 设aii 0;i 1,2,n的SOR 方法迭代公式为x(k1)Lx(k)f.8其中L(D L)1(1)D U,f(D L)1b.矩阵L称为 SOR 方法的迭代矩阵。四、算法程序四、算法程序一Jacobi 迭代法function x,n=jacobi(A,b,x0,eps,varargin)if nargin=3 eps=1.0e-6;M =50000;elseif nargin=eps x0=x;x=B*x0+f;n=n+1;if(n=M)disp(Warning:迭代次数太多，可能不收敛！);return;endend(二)SOR 方法function x,n=SOR(A,b,x0,w,eps,M)if nargin=4 eps=1.0e-6;12/37优质文本 M =50000;elseif nargin4 error returnelseif nargin=5 M =50000;endif(w=2)error;return;endD=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);B=inv(D-L*w)*(1-w)*D+w*U);f=w*inv(D-L*w)*b;x=B*x0+f;n=1;while norm(x-x0)=eps x0=x;13/37优质文本 x=B*x0+f;n=n+1;if(n=M)disp(Warning:迭代次数太多，可能不收敛！);return;endend五、求解过程五、求解过程一第1题1、MATLAB 计算过程在 Command Window 中输入 H=hilb(6);xn=ones(6,1);b=H*xn;x0=zeros(6,1)x0=0 0 014/37优质文本 0 0 0 x1,m1=jacobi(H,b,x0)x1=Inf Inf NaN NaN NaN NaNm1=487 xx11,m11=SOR(H,b,x0,1)xx1=15/37优质文本 0.9999 1.0009 0.9981 0.9974 1.0089 0.9947m11=17406 xx21,m21=SOR(H,b,x0,1.25)xx21=1.0000 1.0003 1.0010 0.9922 1.0126 0.993916/37优质文本m21=16290 xx31,m31=SOR(H,b,x0,1.5)xx31=1.0000 0.9995 1.0052 0.9829 1.0216 0.9907m31=167692、结果分析：从 MATLAB 输出的结果可以看出：当 n=6 时，取初始向量x0=(0,0,0,0,0,0),用 Jacobi 迭代法求解该矩阵方程经过m1=487 次迭代所得结果 x1=(Inf,Inf,NaN,NaN,NaN,NaN),可见17/37优质文本用Jacobi迭代法在这种情况下求解该矩阵方程是不收敛的。用SOR方法求解该矩阵方程，当松弛因子 w=1 时，经过 m11=17406 次迭代所得结果 xx11=(0.9999，1.0009，0.9981，0.9974，1.0089，0.9947),此 结 果 及 精 确 解x*(1,1,1,1,1,1)比 拟 可 得xx11 x*0.0109，可见此结果已经很接近精确解了；当松弛因子w=1.25时，经 过m21=16290次 迭 代 所 得 结 果xx21=(1.0000,1.0003,1.0010,0.9922,1.0126,0.9939),此结果及精确解x*(1,1,1,1,1,1)比拟可得xx21 x*0.0160，可见此结果没有w=1 所得的解接近精确解，但迭代次数比它少；当松弛因子w=1.5时，经过 m31=16769 次迭代所得结果 xx31=(1.0000，0.9995，1.0052，0.9829，1.0216，0.9907),此结果及精确解x*(1,1,1,1,1,1)比拟可得xx31 x*0.0296，可见此结果没有 w=1 和 w=1.25 所得的解接近精确解，迭代次数虽然比 w=1 时少，但比 w=1.25 多。二第2题1、MATLAB 计算过程(1)取 n=8在 Command Window 中输入 H=hilb(8);xn=ones(8,1);b=H*xn;18/37优质文本 x0=zeros(8,1)x0=0 0 0 0 0 0 0 0 x2,m2=jacobi(H,b,x0)x2=-Inf NaN NaN NaN NaN19/37优质文本 NaN NaN NaNm2=396 xx12,m12=SOR(H,b,x0,1)xx12=1.0001 0.9974 1.0136 0.9794 0.9982 1.0141 1.0101 0.9870m12=20/37优质文本 8342 xx22,m22=SOR(H,b,x0,1.25)xx22=1.0002 0.9964 1.0179 0.9726 1.0031 1.0130 1.0091 0.9877m22=17436 xx32,m32=SOR(H,b,x0,1.5)xx32=21/37优质文本 1.0001 0.9976 1.0119 0.9812 1.0043 1.0055 1.0079 0.9915m32=37921(2)取 n=10在 Command Window 中输入 H=hilb(10);x1=ones(10,1);b=H*x1;x0=zeros(10,1)x0=22/37优质文本 0 0 0 0 0 0 0 0 0 0 x3,m3=jacobi(H,b,x0)x3=Inf Inf Inf NaN NaN NaN23/37优质文本 NaN NaN NaN NaNm3=347 xx13,m13=SOR(H,b,x0,1)xx13=1.0001 0.9982 1.0056 0.9993 0.9916 0.9968 1.0052 1.0087 1.003924/37优质文本 0.9904m13=26951 xx23,m23=SOR(H,b,x0,1.25)xx23=1.0001 0.9989 1.0022 1.0042 0.9903 0.9965 1.0042 1.0080 1.0039 0.9916m23=25/37优质文本 30092 xx33,m33=SOR(H,b,x0,1.5)xx33=1.0000 0.9997 0.9986 1.0109 0.9850 0.9993 1.0020 1.0086 1.0033 0.9924m33=321742、结果分析：从MATLAB输出的结果可以看出：(1)用Jacobi迭代法求解该矩阵方程当n=8时，取初始向量x0=(0,0,0,0,0,0,0,0),经过m2=396次26/37优质文本迭代所得结果x2=(-Inf，NaN，NaN，NaN，NaN，NaN，NaN，NaN)；当n=10时，取初始向量x0=(0,0,0,0,0,0,0,0,0,0),经过m3=347次迭代所得结果x3=(Inf，Inf，Inf，NaN，NaN，NaN，NaN，NaN，NaN，NaN)，比照1题中所得结果可见无论n确何值用Jacobi迭代法在这种情况下求解该矩阵方程是不收敛的，其原因为其迭代矩阵B D1(L U)(其中D=diag(diag(H),L=-tril(H,-1),U=-triu(H,1)的谱半径2用SOR方法求解该矩阵方程：当(B)4.3085不满足(B)1。n=8时，取初始向量x0=(0,0,0,0,0,0,0,0)当松弛因子w=1时，经过m12=8342次迭代所得结果xx12=(1.0001，0.9974，1.0136，0.9794，0.9982，1.0141，1.0101，0.9870),此结果及精确解x*(1,1,1,1,1,1,1,1)比拟可得xx12 x*0.0330；当松弛因子w=1.25时，经过m22=17436次迭代所得结果xx22=(1.0002，0.9964，1.0179，0.9726，1.0031，1.0130，1.0091，0.9877),此结果及精确解x*(1,1,1,1,1,1,1,1)比拟可得xx22 x*0.0387，；当松弛因子w=1.5时，经过m32=16769次迭代所得结果xx32=(1.0001，0.9976，1.0119，0.9812，1.0043，1.0055，1.0079，0.9915),此结果及精确解x*(1,1,1,1,1,1,1,1)比拟可得xx32 x*0.0262；当n=10时，取初始向量x0=(0,0,0,0,0,0,0,0,0,0)当松弛因子w=1时，经过m13=2695次迭代所得结果xx13=(1.0001，0.9982，1.0056，0.9993，0.9916，27/37优质文本0.9968，1.0052，1.0087，1.0039，0.9904),此结果及精确解x*(1,1,1,1,1,1,1,1,1,1)比拟可得xx13 x*0.0181；当松弛因子w=1.25时，经过m23=30092次迭代所得结果xx23=(1.0001，0.9989，1.0022，1.0042，0.9903，0.9965，1.0042，1.0080，1.0039，0.9916),此结果及精确解x*(1,1,1,1,1,1,1,1,1,1)比拟可得xx23 x*0.0173，；当松弛因子w=1.5时，经过m33=32174次迭代所得结果xx33=(1.0000，0.9997，0.9986，1.0109，0.9850，0.9993，1.0020，1.0086，1.0033，0.9924),此结果及精确解x*(1,1,1,1,1,1,1,1,1,1)比拟可得xx33 x*0.0222,比照n=8,n=10时，w取相同值时的结果，当n=10时的迭代步骤比n=8都少，且所得解比起n=8也精确。三第3题在1、2题中所得的结果看来无论 n 确何值用 Jacobi 迭代法在这种情况下求解该矩阵方程是不收敛的,而用 SOR 方法解题时，虽然迭代收敛，但需要迭代很屡次，针对这次问题，以下尝试采用另外一种古典迭代法 Gauss-Seidel 法解决此题：1、Gauss-Seidel 法程序function x,n=gauseidel(A,b,x0,eps,M)if nargin=3 eps=1.0e-6;M =50000;28/37优质文本elseif nargin=4 M =50000;elseif nargin=eps x0=x;x=G*x0+f;n=n+1;if(n=M)disp(Warning:迭代次数太多，可能不收敛！);return;29/37优质文本 endend2、MATLAB 计算过程在 Command Window 中输入 H=hilb(6);x1=ones(6,1);b=H*x1;x0=zeros(6,1)x0=0 0 0 0 0 0 xxx,mnn=gauseidel(H,b,x0)xxx=0.9999 1.0009 0.998130/37优质文本 0.9974 1.0089 0.9947mnn=17406由此题结果可见，用此方法解题，及 SOR 方法，当松弛因子w=1 时所得结果相同，正好验证了上面所说的：当松弛因子 w=1 时，SOR 方法就是 Gauss-Seidel 法。虽然此方法并不比 SOR 方法更简便，但对于不收敛的 Jacobi 迭代法来说，此方法是可以解出此题的。六、心得体会六、心得体会通过这次专业综合训练，对题目的研究求解，对 Hilbert 矩阵方程组及其特点有了一定的了解，熟练掌握运用了古典迭代法中的Jacobi 迭代法、SOR 方法，并对 Gauss-Seidel 法也做了一定的了解，对这些这些古典迭代法的特性有了比拟系统的认识。对 MATLAB软件的运用也更熟练了。七、参考文献七、参考文献1谢进，李大美.MATLAB 及计算方法实验M.武汉大学出版社，2016.2周品，何正风.MATLAB 数值分析M.机械工业出版社，2016.131/37优质文本3林雪松，周婧，林德新.MATLAB7.0 应用集锦M.机械工业出版社，20064(美)John H.Mathews,Kurtis D.Frink 数值方法M.电子工业出版社，20055李庆扬，王能超，易大义.数值分析M.华中科技大学出版社，2006.132/37
