数值线性代数课设资料

1、数值线性代数课程设计报告姓名:陶英学号:081410124任课教师:杨熙南京航空航天大学2016年6月22日求解线性方程组的三种迭代法及其结果比较摘要当今的环境下，数值计算越来越依赖于计算机。大规模科学计算和工程技术 中许多问题的解决，最终归结为大型稀疏线性方程组的求解，其求解时间在整个问题求解时间中占有很大的比重，有的甚至达到 80%。由于现今科学研究和大型 项目中各种复杂的可以对计算精度和计算速度的要求越来越高。因此，作为大规模科学计算基础的线性代数方程组的高效数值求解引起了人们的普遍关注。这种方程组的求解一般采用迭代法。关于迭代法，是有很多种解决公式的： Jacobi，G-S和超松弛迭代

2、法。这三 种方法的原理大致相同，Jacobi需要给定初向量，G-S则需要给定初值，超松弛 法是对Guass-Seidel迭代法的加权平均改造。而本文则是对大型稀疏线性方程组 迭代求解与三种迭代法（Jacobi, Gauss-Seide和超松弛迭代法）的收敛速度与精 确解的误差比较做出研究。关键词：Jacobi迭代法；Gauss-Seidel迭代法；SOR迭代法；线性方程组1方法与理论的叙述1.1迭代法简介I. Jacobi迭代法：对于非奇异线性方程组 Ax=b,令A=D-L-U,其中(2.2)D =(iiag(airan g).z =0柑 2L03-為i-如0一如 一- an.n-lu =&#

3、176; -au -0 -如0 一 一0 _则原方程组可改写为：xBjX + gj 其中Bj = D(L + U).gj = Dlb,给定初始向量：由(2.2 )可以构造迭代公式：B了兀卜+g,k = 1*2*3 *其分量形式为：2. Guass-Seidel 迭代法：类似于Jacobi迭代法，给定初值:令% =（!>£尸5 g俯-工尸打则得到Guass-Seidel公式：其分量形式为:3. 超松弛迭代法（SOR迭代法）：SOF迭代法是对Guass-Seidel迭代法的加权平均改造，即工=(1-少)严乜+获¥ 为Guass-Seidel迭代解，即严)=(1 一+ 口(

4、矿仏")+ 矿1叭它的分量形式为：1 j-JJI尹=(1+少丄(二。眄+二。内+处aa j=ij=»-i其中3称为松弛因子，当3 >1时称为超松弛；当3 <1时叫低松弛；3 =1时就是Guass-Seidel 迭代。上述三种经典迭代法收敛的充分必要条件是迭代矩阵谱半径小于1。谱半径不易求解，而在一定条件下，通过系数矩阵A的性质可判断迭代法的收敛性。定理1:若系数矩阵A是严格对角占优或不可约对角占优，则Jacobi迭代法和Gauss-Seide迭代法均收敛。定理2：(1) SOR迭代法收敛的必要条件是 0<w<2;(2) 若系数矩阵A严格对角占优或不可

5、约对角占优且 0<w<1,则SOR迭代法收敛。w=1时，SOR迭代法退化为Gauss-Seide迭代法2数值结果2.1问题考虑两点边值问题：了- d2y dy2a,0 ： a : 1“ dx2 dx、y（0）=0,y（1）=1.1x容易知道它的精确解为：y二：/ .（1 -e：） ax为了将微分方程离散，把0,1区间n等分，令h=1/n，x =ih,i =1,2,n-1，得到差分方程（； h）yi 1 - （2 ； h）% y二二ah2，从而得到迭代方程组的系数矩阵A对；=1, a=1/2,n=100,分别用jacobi，G-S,超松弛迭代法分别求线性方程组的解，要求4位有效数字，

6、然后比较与精确解的误差。对；=0.1，； =0.01，； =0.001，考虑同样问题。1. 方程的表示及存储由于本题中线性方程组的系数矩阵为三对角矩阵，所以可以采用紧缩方法存储,0-（2拦+方）£A-hs -(le + h)c-h£- (2"方)£一（2左+月）0然后在矩阵乘法时对下标处理一下即可。但是考虑到三种迭代方法的一般性，且本题中n=200并不是很大，所以实验中并没有采用紧缩存储， 而是采用了 直接存储。2. 边值条件的处理由于差分得到的方程组的第一行和最后一行中分别出现了边值y（0）与y（1）作为常数项，因此要在常向量的第一项和最后一项作一些修

7、改：h(l) = an y(0) * s-1) = ah' -+3. 迭代终止条件首先确定要求的精度tol ，我们希望当忑爲Xrol则停止迭代。对于迭代格式'-"，若7 _ 一且f d，则迭代序列的II第k次近似解和精确解之间有估计式。由题目要求知我们需要有，而由上面的迭代估计，只要 |- Xi IK io" +时，SOR方法的迭代速度最快。 Matlab 实现： 迭代矩阵是n-1阶的，不是n阶； 等号右端向量b的最后一项，不是ahA2，而是ahA2-eps-h I %】一耳 IH1° 41_9，即q 即可。而本题中q可近似取为-，因此最后令迭代终

8、止条件为q4.S0R迭代中最佳松弛因子的选取由于SOR迭代法的效果和其松弛因子 w的选取有关，所以有必要选取合适的松弛因子。当选择最佳松弛因子2.2精确解-e ；) ax带入 a=1/2, : =1代码：>> clear>> x=li nspace(0,1);plot(x,truy,'g','Li neWidth',1.5);hold on;Grid图:ure 12.3三种迭代法Jacobi法：代码见附录Eps=1结果：迭代次数k: 22273结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）0.8268 (L 8354 0.8

9、UC 0.8525A. 86090.86940.877SQ.88S20.89450.9028o.$ino.diasQu 92750.93570.94380.95200. AMOoiseeia.sTei0.98410.992122273Eps=0.1结果：迭代次数k: 8753结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）<L BS97CL 89486 B9980.90480.90980.91490.91996 92490.9299(L 93490.93990.9449Q« M990.9550仇 96000.9650<L &7000.97500.980

10、00.9B500.9900山 99508753Eps=0.01结果：迭代次数k: 661结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）0,89500+ 90000.90500.91000,91500.9200Q.9250Q.93000.93500.94000.94500.95000.95500. 96000.96500l 97000- 97500.98000.98500+990Q0.99 禎G-S迭代法：代码见附录Eps=1结果：迭代次数k: 11125结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）low HeliCurrent FCommanc60.Q0.CL0.

11、0JQ060.心 New to oJCL0+0.97610.9841(L 992111125Eps=0.1结果：迭代次数k: 4394结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）0, £99?Q. 904B0.90950. SUS0.91»90. 924»Q. 9E990. 9349Q. 93990.944B0. 9499O,$55O0.9600Q. 96 BO0. 97000. 975&0.9SO00. 9B5Q0. S9000. 99504394Eps=O.O1结果：迭代次数k: 379结果与精确解的比较图（绿色粗线是精确解，黑色细线是

12、迭代结果）0. 90000.905Q0.91000.91500. 9200<92500.93000.93500. 94000.34500,95000.95500.96000.96500.97000.97500. 98000.9E500.95000. 9950超松弛法：代码见附录Eps=1 w=1.56结果：迭代次数k: 3503结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）LL 451 Q：0,82630. 83640.84400. 8525 £6090, 86940-87786 88620, 89450. 902 B0,91110.91930 92750,935

13、70.94380+ 95200, 96000. S6B10,97610. 9841CL 9921£ Figrure 2Ele Edit 如申 Insert Jools desktop Jindcw fc|ulp d d| S O ® '-G-层J 口 目旦3503Eps=0.1 w=1.56结果：迭代次数k: 1369结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）0,88870.S9960.904S0.90880.91490.91990.92490.92990.$34$0.93QS0,94490.94990.95500.9600Ou 9660Q.970

14、00.97500.98000.93500.991300.9950Eps=0.01 w=1.56 结果： 迭代次数k: 131结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）0,90000.90500. 91000. 91600. 92000, S2500. 9300(L 93500. 94000.94500. 95000. 9550 96000. 96500. 97000. 97500. 98000. 9S500. 99000. 9960 Figure 2File £dit 里 iew Insert Iools JQektop JMindw Htlp己彳k毀、凸動县|口用

15、旦3分析讨论及心得体会3.1三种方法的比较Jacobi、G-S、超松弛法，三者都能够取得对精确解的良好逼近，但是，在相同 的精度条件下，三者的收敛速度是不一样的，jacobivG-Sv超松弛，也就是说， 在迭代次数相同的条件下，精度：jacobivG-Sv超松弛。3.2心得体会这次课程设计，平时感觉挺简单的那些枯燥单调的代码和数学公式，真正到了自己运用的时候却无从下手，但是，解决问题的过程恰是不断学习的过程：数学算法转换为代码的过程要对题目有深入的了解，然后对程序函数定义还要有一定的掌握能力，通过这个的过程让我巩固了自己的数学知识，对数学专业知识和MATLAB勺操作有了更深的体会。课程设计中遇

16、到的问题只凭自己苦思冥想是不能全部解决的， 这是同学老师的 建议和网络给了我很大的帮助。遇到自己解决不了的问题时，多多向老师同学请 教，或许问题就能迎刃而解。4参考文献1 徐树方数值线性代数.北京:北京大学出版社,1995.2 马昌凤现代数值分析.北京:国防工业出版社.2013.3 刘春凤，米翠兰实用数值分析教程北京冶金工业出版社.20065附录源代码I. Jacobi :fun cti on y,k=jacobi2(a,eps,h,delta) n=1.0/h;A=on es( n-1);y=zeros( n-1,1);z=zeros( n-1,1);k=0;for i=1: n-1for

17、j=1: n-1 A(i,j)=0; end end for i=1: n-1 A(i,i)=-(2*eps+h);endfor i=1: n-1for j=1: n-1 if i=j+1A(i,j)=eps;end if i=j-1A(i,j)=eps+h;endendendb=zeros( n-1,1);for i=1: n-2b(i,1)=a*hA2;endb(n-1,1)=a*hA2-eps-h;D=zeros( n-1);for i=1: n-1D(i,i)=A(i,i);endL=zeros( n-1);for i=1: n-1for j=1: n-1if i>jL(i,j)

18、=-A(i,j);endendendU=zeros( n-1);for i=1: n-1for j=1: n-1if i<jU(i,j)=-A(i,j);endendendB=D(L+U);g=Db;while 1z=B*y+g;if norm(z-y,inf)<delta break;endy=z;k=k+1;endx=li nspace(0,1);truy=(1-a)/(1-exp(-1/eps)*(1-exp(-x./eps)+x.*a; figure;plot(100*x,truy, 'g' , 'LineWidth' ,5);hold on

19、;gridhold on;plot(y, 'b')2. G-S: fun cti on y,k=gs2(a,eps,h,delta) n=1.0/h;A=on es( n-1);y=zeros( n-1,1);z=zeros( n-1,1);k=0;for i=1: n-1for j=1: n-1 A(i,j)=0; end end for i=1: n-1 A(i,i)=-(2*eps+h);endfor i=1: n-1for j=1: n-1 if i=j+1 A(i,j)=eps;end if i=j-1A(i,j)=eps+h;endendendb=zeros( n-

20、1,1);for i=1: n-2b(i,1)=a*hA2;endb(n-1,1)=a*hA2-eps-h;D=zeros( n-1);for i=1: n-1D(i,i)=A(i,i);endL=zeros( n-1);for i=1: n-1for j=1: n-1if i>j L(i,j)=-A(i,j);endendendU=zeros( n-1);for i=1: n-1for j=1: n-1if ivjU(i,j)=-A(i,j);endendendB=D(L+U);g=Db;while 1z=(D-L)U*y+(D-L)b;if norm(z-y,inf)<delta break;endy=z;k=k+1;endx=li nspace(0,1);truy=(1-a)/(1-exp(-1/eps)*(1-exp(-x./eps)+x.*a; figure;plot(100*x,truy, 'g' , 'LineWidth' ,5);hold on;gridhold on;plot(y, 'b')3.SORfun cti on y,k=sor(a,eps,h,delta,w) n=1.0/h;A=on es( n-1);y=zeros( n-1,1);z=zeros