古典迭代法中的Jacobi迭代法和SOR迭代法求解Hilbert矩阵方程组

1、燕山大学课程设计说明书题目：直接法求解Hilbert 矩阵方程组学院（系）：理学院年级专业：学号：学生姓名：指导教师：教师职称：教授燕山大学课程设计（论文）任务书院（系）：理学院基层教学单位：学 号学生姓名专业（班级）设计题目古典迭代法求解 Hilbert 矩阵方程组设分别利用古典迭代法中的Jacobi 迭代法和 SOR 迭代法，求解 Hilbert 矩阵方程计组，并与准确解进行比较，得出结论，然后针对不同的n 值重复进行计算并比技较所得的结果。术参数写出相应的 Matlab 程序，执行程序得出的结果。设通过结果对比分析并进行理论分析。计要求四周的时间掌握 Hilbert 矩阵的性质、 Hil

2、bert 矩阵方程组的特点及其相关的求解方法，研究古典迭代法中的Jacobi 迭代法和 SOR 迭代法并用这两种方法求解工Hilbert 矩阵方程组，分析这两种办法对解Hilbert矩阵方程组的效果如何；作熟悉使用 Matlab 软件，上机求解给定题目。量第一周： 查阅资料， 基本了解 Jacobi 迭代法和 SOR 迭代法的基本原理、解题步工骤等，主要对Jacobi 迭代法进行分析；作第二周至第三周： 根据其他组员编写的Jacobi 迭代法的Matlab 程序进行具体题计目的计算并与其他同学负责的部分进行综合分析；划第四周：整理和撰写报告，进行答辩。参1 谢进，李大美 . MATLAB与计算

3、方法实验 M.武汉大学出版社，2009.2 周品，何正风 .MATLAB数值分析 M. 机械工业出版社， 2009.1考3 林雪松，周婧，林德新.MATLAB7.0 应用集锦 M. 机械工业出版社， 2006资4( 美 )John H.Mathews,Kurtis D.Frink 数值方法 M. 电子工业出版社，2005料5 李庆扬，王能超，易大义 .数值分析 M. 华中科技大学出版社， 2006.1指导教师签字基层教学单位主任签字说明：此表一式四份，学生、指导教师、基层教学单位、系部各一份。2012 年 01月 09日燕山大学课程设计评审意见表指导教师评语：成绩：指导教师：年月日答辩小组评语

4、：成绩：评阅人：年月日课程设计总成绩：答辩小组成员签字：年月日燕山大学课程设计说明书目录摘要 2 一、课题背景简介 3 二、数学模型描述 3三、算法原理 4 (一) Hilbert 矩阵的来源及其正定性 .4(二)迭代法的基本思想 . 5 （ 1）Jacobi 迭代法 .6 （ 2）SOR 方法 7四、算法程序 8(一)Jacobi 迭代法 8(二) SOR 方法 9五、求解过程 10（一）第（ 1）题 10（二）第（ 2）题 12（三）第（ 3）题18六、心得体会 20七、参考文献 20共20页第1页燕山大学课程设计说明书摘要论文旨在阐述求解 Hilbert 矩阵方程方法中的古典迭代法应用于

5、求解特殊矩阵的过程 ,其中主要包括 Jacobi迭代法和 SOR 迭代法。文中从方法原理、 算法的基本应用到理论的归纳与扩展。通过这一例题，可进一步领会数值分析的实际应用。同时为了考察读者对古典迭代法的学习灵活度，例题的解答需要用到数值方法在MATLAB 中的实现，包括数值方法在MATLAB 中的函数实现， MATLAB编程和 MATLAB绘图，使读者在上机练习中加深对数值分析算法原理的理解，通过对算法的思考和理论分析，掌握MATLAB 的使用。关键词 ：Hilbert 矩阵古典迭代法Jacobi 迭代法 SOR 迭代法共20页第2页燕山大学课程设计说明书一、课题背景简介随着电子计算机的出现和

6、迅速发展，在各门自然科学和工程技术科学的发展中，科学计算已经成为平行于理论分析和科学实验的第三种科学手段。数值计算是科学计算中的一个不可缺少的环节，而在数值计算中，一类很重要的问题就是线性方程组的求解。另外，数学和物理以及力学等学科和工程技术中许多问题的最终解决都归结为解一个或一些大型系数矩阵的线性方程组。所以，大型线性方程组的求解是大规模科学与计算的核心，许多作者都对此作了研究。对线性方程组 Ax b 的求解，主要有直接求解法和迭代法求解。对于阶数不太高的线性方程组，用直接法比较方便，高斯消元法和克兰姆法则是直接解法里面最重要的解法。但是，随着科学技术的飞速发展，需要求解的问题的规模越来越大

7、，迭代法已取代直接法成为求解大型线性组的最重要的一类方法。此时，迭代格式的收敛性和收敛速度成为一个很重要的问题，成为人们关注的焦点。不收敛的格式当然不能用，虽然收敛但是收敛的很慢的格式，不仅是人工和机器的时间比较浪费，而且还不一定能解出结果，实际应用价值太小。 通常用的迭代法有 Jacobi，Gauss-Seidel 等古典迭代法，还有SOR， AOR 以及 SAOR， SSOR 等迭代法。这些迭代法的提出对大型线性方程组的求解提供了一条快速有效的途径。但是这些迭代法相对来说使用条件很苛刻，或者很繁琐。此外当迭代矩阵的谱半径比较大，尤其接近 1 时，迭代速度将会很慢，极大的影响了方程组求解的时

8、效性，从而使它们的应用受到了一定的限制。MATLAB 是矩阵实验室 （Matrix Laboratory ）的简称，是美国 MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括 MATLAB 和 Simulink 两大部分。MATLAB 和 Mathematica、Maple 并称为三大数学软件。 它在数学类科技应用软件中在数值计算方面首屈一指。 MATLAB 可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设

9、计与分析等领域。MATLAB 的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB 来解算问题要比用C，FORTRAN 等语言完成相同的事情简捷得多，并且mathwork 也吸收了像 Maple 等软件的优点 ,使MATLAB 成为一个强大的数学软件。 在新的版本中也加入了对C，FORTRAN ，C+ ，JAVA 的支持。可以直接调用,用户也可以将自己编写的实用程序导入到 MATLAB 函数库中方便自己以后调用，此外许多的 MATLAB 爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。MATLAB 的优势：共20页第3页燕山大学课程设计说明书（1）

10、友好的工作平台和编程环境； （ 2）简单易用的程序语言；（3）强大的科学计算机数据处理能力； （4）出色的图形处理功能； （ 5）应用广泛的模块集合工具箱；（6）实用的程序接口和发布平台； （ 7）应用软件开发（包括用户界面）二、数学模型描述：古典的迭代法求解Hilbert 矩阵方程组考虑方程组 Hxb ,其中系数矩阵为：H ( hij ) R n n , hiji1,i , j 1,2, , nj 1假设由准确解 x(1,1,1)TRn 确定向量 b 。（1）选择 n=6，分别用 Jacobi 迭代法和 SOR 迭代法（1,1.25,1.5 等）求解。将计算结果与准确解进行比较。（2）逐步增

11、大 n(n=8,10,),重复（ 1）的计算。比较结果，结果说明什么？（3）对于计算过程中发现的问题，尝试用其他方法加以解决。三、算法原理 :(一) Hilbert 矩阵的来源及其正定性设 f ( x) C0,1 ,求 n 次多项式P( x) a0a1 x a2 x 2an x n使得L12 dx P( x) f (x)0取最小值。其中， P( x) 称为函数 f ( x) 的最佳平方逼近多项式。由极值必要条件，对1na j x jf ( x) 2 dxL (a0 , a1, , an )0j 0中的变量求导数，得L1na j x jn1xi j12 xif (x)dx2aj0dx 2 xi

12、f ( x) dxai0j0j00令其为零，得方程组共20页第4页燕山大学课程设计说明书n11a j(i0,1, n)xi f (x)dxj 0 ij10令 bi1xi f (x)dx ，将方程组写为矩阵形式011/ 21/( n1)a0b01/ 21/ 31/(n2)a1b11/(n1)1/( n2)1/(2n1)anbn方程组的系数矩阵就是著名的Hilbert 矩阵。多项式可表示为内积形式1xP( x) a0a1an xn所以1a0 P( x) 2 a0 a1x1xxn a1an x nan积分，得11/ 21/(n1)a012 dx a0 a1an 1/ 21/ 31/( n 2)a1

13、P( x)01/( n 1)1/(n2)1/( 2n1)an由此可知， Hilbert 矩阵是对称正定矩阵。另外， Hilbert 矩阵是著名的病态矩阵。(二)迭代法的基本思想对于大型方程组Ax b（1）其中n nnTAR ， ，bR 。假设 A 非奇异，方程组有唯一解 x( x , x, x) 。x12n将 Axb 变形为等价的方程组共20页第5页燕山大学课程设计说明书xBxf（2）由此建立迭代公式x (k 1)Bx( k )f , k0,1,2,（ 3）给定初始向量 x (0) ，按此公式计算得近似解向量序列x( k ) ，称此方法为迭代法。若对任意 x( 0) ，当迭代次数无限增加时，序

14、列x( k) 都有相同的极限 x ，即 lim x (k )x ，也就是说 lim x( k)x ， i1,2, n 。显然有kkxBxf（4）这是说迭代公式是收敛的，否则说发散的，称迭代格式（3）中的矩阵 B 为迭代矩阵，对于不同的迭代矩阵得到不同的迭代格式。（1）Jacobi 迭代法设有方程组naij xjbi( i1,2, n ),j1记作Axb ,（1）A 为非奇 异 矩阵 aij0(i1,2, n)。 将 A分解为A D LU ,其中110aa22a210DLa31a320ann，an1an 2an ,n 10，0a12a13a1n0a23a2nU0an1,n0共20页第6页燕山大学

15、课程设计说明书将式（ 1）的第 i(i1,2, n) 个方程用 aij 去除再移项，得到等价方程组1n, n) ，（2）xi(biaij x j )(i1,2,aiij 1j i简记作xB0 xf ,其中 B0I D 1 A D 1 ( L U ) , f D 1b 。对方程组（ 2）应用迭代法，得到解式（1）的 Jacobi 迭代公式x( 0)0( 0)(o) T(初始向量 ，( x1, x2, , xn )( k 1)1n( k)（3），xi(biaij x j )aiii 1j i其中 x( k)(x1k , x2( k) , , xn( k) )T 为第 k 次迭代向量。设 x ( k

16、 ) 已经算出，由式 (3)可计算下一次迭代向量 x( k 1)（ k 0,1,2,; i 1,2, n ）。x(0 )(初始向量 ),显然迭代公式（ 3）的矩阵形式为B0 x( k)x( k 1)f ,其中 B0 称为 Jacobi 迭代矩阵。（2）SOR 方法设有方程组Axb,（ 1）其中 AR nn 为非奇异矩阵，且设 aii0（ i1,2,n ），分解 A 为A DL U 。（ 2）设已知第 k次迭 代向 量 x (k )， 及第 k+1次迭 代 向量 x( k 1)的分 量 x(j k 1 )（ j 1,2, , i1）,要求计算分量 xi( k 1) 。首先用 Gauss-Seid

17、el迭代法定义辅助量(k 1)1i 1n(k 1)( k )(i1,2, n) ，（3）xi(biaij x jaij x jaiij 1j i 1共20页第7页燕山大学课程设计说明书( k 1 )( k )( k 1)某个平均值（及加权平均） ，即再把 xi取为 xi与 xi(k !)(1(k ) ( k ! )( k) (k !)( k )（4）xi)xixixi(xixi用式 (3)代入式（ 4）即得到解方程组Axb 的逐次超松弛迭代公式i1nxi( k !)xi( k)aii(biaij x (jk !)aij x(j k) ),（ 5）j 1jix( k )(x1(k ) , x2(

18、k ) , xn(k ) )T (k0,1,; i1,2, n),其中称为松弛因子，或写为xi(k!)xi(k )xi(k0,1,;i1,2, n),i1aij x(j k 1)naij x(j k ) ).（ 6）xiaii(bij1ji显然，当1时，解式（1）的 SOR 方法就是 Gauss-Seidel迭代法。当1时，称式（ 5）为低松弛法，当1时称式（ 5）为超松弛法。SOR 迭代公式的矩阵形式，迭代公式（5）也可写为aii xi(k1)xi(k)i1x(j k1)naij x(j k )(1)aii(biaij)(i1,2, n) （7）j1j i1用式 ADLU ，则Dx (k1)

19、(bLx (k1)Ux( k) )(1) Dx (k )，即DL x(k 1)(1)DU)x( k)b.()显然对于任何一个值， (DL) 非奇异（由设 aii0; i1,2, n ），于是x (k !)(DL) 1 (1)DU x (k )( DL) 1 b.也就是说，解式（ 1）（设 aii0;i1,2,n ）的 SOR 方法迭代公式为x (k 1)L x (k )f .（8）其中L( DL) 1(1)DU , f( DL) 1b .矩阵 L 称为 SOR 方法的迭代矩阵。共20页第8页燕山大学课程设计说明书四、算法程序（一） Jacobi 迭代法function x,n=jacobi(A

20、,b,x0,eps,varargin)if nargin=3eps= 1.0e-6;M= 50000;elseif nargin<3errorreturnelseif nargin =5M= varargin1;endD=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);B=D(L+U);f=Db;x=B*x0+f;n=1;while norm(x-x0)>=epsx0=x;x=B*x0+f;n=n+1;if (n>=M)disp('Warning: 迭代次数太多，可能不收敛！ '); return;endend(二) SOR 方

21、法function x,n=SOR(A,b,x0,w,eps,M)if nargin=4eps= 1.0e-6;M= 50000;elseif nargin<4errorreturnelseif nargin =5M= 50000;共20页第9页燕山大学课程设计说明书endif(w<=0 | 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)>=epsx

22、0=x;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 =000000共20页第10页燕山大学课程设计说明书>> x1,m1=jacobi(H,b,x0)x1 =InfInfNaNNaNNaNNaNm1 =4

23、87>> xx11,m11=SOR(H,b,x0,1) xx1 =0.99991.00090.99810.99741.00890.9947m11 =17406>> xx21,m21=SOR(H,b,x0,1.25) xx21 =1.00001.0003共20页第11页燕山大学课程设计说明书1.00100.99221.01260.9939m21 =16290>> xx31,m31=SOR(H,b,x0,1.5) xx31 =1.00000.99951.00520.98291.02160.9907m31 =167692、结果分析：从 MATLAB 输出的结果可以

24、看出：当 n=6 时，取初始向量 x0=(0,0,0,0,0,0),用 Jacobi 迭 代法 求解 该矩 阵方 程经 过 m1=487 次 迭 代 所 得结 果 x1=(Inf,Inf,NaN,NaN,NaN,NaN) ,可见用 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)' 比较可得xx11x*0.0109 ，可见此结果已经很接近精确

25、解了；当松弛因子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)' 比较可得xx21x*0.0160 ，可见此结果没有w=1 所得的解接近精确解，但迭代次数比它少；当松弛因子w=1.5 时，经过 m31=16769 次共20页第12页燕山大学课程设计说明书迭代所得结果 xx31 =( 1.0000 ，0.9995，1.0052，0.9829，1.0216，0.9907),此结果与精确解 x* (1,1,1,1,1,1)'

26、比较可得 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;>> x0=zeros(8,1)x0 =00000000>> x2,m2=jacobi(H,b,x0)x2 =-InfNaNNaNNaNNaNNaNNaNNaN共20页第13页燕山大学课程设计说明书m2

27、 =396>> xx12,m12=SOR(H,b,x0,1) xx12 =1.00010.99741.01360.97940.99821.01411.01010.9870m12 =8342>> xx22,m22=SOR(H,b,x0,1.25) xx22 =1.00020.99641.01790.97261.00311.01301.00910.9877m22 =17436共20页第14页燕山大学课程设计说明书>> xx32,m32=SOR(H,b,x0,1.5) xx32 =1.00010.99761.01190.98121.00431.00551.0079

28、0.9915m32 =37921(2)取 n=10在 Command Window 中输入>> H=hilb(10);>> x1=ones(10,1);>> b=H*x1;>> x0=zeros(10,1)x0 =0000000000>> x3,m3=jacobi(H,b,x0)共20页第15页燕山大学课程设计说明书x3 =InfInfInfNaNNaNNaNNaNNaNNaNNaNm3 =347>> xx13,m13=SOR(H,b,x0,1)xx13 =1.00010.99821.00560.99930.99160.9

29、9681.00521.00871.00390.9904m13 =26951共20页第16页燕山大学课程设计说明书>> xx23,m23=SOR(H,b,x0,1.25) xx23 =1.00010.99891.00221.00420.99030.99651.00421.00801.00390.9916m23 =30092>> xx33,m33=SOR(H,b,x0,1.5) xx33 =1.00000.99970.99861.01090.98500.99931.00201.00861.00330.9924m33 =321742、结果分析：从MATLAB 输出的结果可以看

30、出： (1)用Jacobi迭代法求解该矩阵方程当 n=8 时，取初始向量 x0=(0,0,0,0,0,0,0,0),经过 m2=396次迭代所得结果 x2 =( -Inf ，NaN，NaN，NaN，NaN， NaN， NaN，NaN)； 当n=10时，取初始向量共20页第17页燕山大学课程设计说明书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迭代法在这种情况下求解该矩阵方程是不收敛的，其原因为

31、其迭代矩阵 BD 1 ( LU ) (其中 D= diag(diag(H),L=-tril(H,-1),U=-triu(H,1)的谱半径(B)4.3085 不满足(B)1。（2）用SOR方法求解该矩阵方程：当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)' 比较可得xx12x*0.0330 ；当松弛因子 w=1.25时，经过 m22=

32、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)' 比较可得xx22x*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)' 比较可得xx32x*0.0262 ；当 n=10时，取初始向量 x0=(

33、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，0.9968，1.0052，1.0087，1.0039，0.9904),此结果与精确解x*(1,1,1,1,1,1,1,1,1,1)' 比较可得xx13x*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),此

34、结果与精确解 x* (1,1,1,1,1,1,1,1,1,1)'比较可得xx23x*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)' 比较可得xx33x*0.0222 ,对比 n=8,n=10时， w取相同值时的结果，当n=10时的迭代步骤比 n=8都少，且所得解比起 n=8也精确。（三）第（ 3）题在（ 1）、（ 2）题中所得的结果看来无论 n 确何值用 Jacobi 迭代法在这种情况下求解该矩阵方程是不收敛的 ,而用 SOR 方法解题时，虽然迭代收敛， 但需要迭代很多次， 针对这次问题， 以下尝试采用另外