迭代法试验报告_米瑞琪

1、数值实验报告n实验名称迭代法求解Poisson方程实验时间2016年4月10日姓名米瑞琪班级信息1303学心1309010304成绩一、实验目的，内容1、理解并掌握三类迭代方法(Jacobi, Gauss-Siedel,SOR)迭代方法的构造原理：2、了解矩阵在matlab中不同的存储方式会导致不同的计算效率，学会采用最高效的方式存储矩阵；3、学会在计算机上实现迭代法，并比较不同方法的效率与误差;二、算法描述(一)Jacobi 迭代对于线性方程组：Au = b (1)为了构造迭代格式，可将上式改写为：u = Tu + d (2)假定A的对角元均不为0,将A分裂为：A = D-B (3)其中D为

2、：D = diag(Aii，4zz，433,，4”)(4)则(1)式可写为：Du = Bu + b (5)形成Jacobi迭代：uk+1 =+ D-lb (6)为了在迭代格式中不出现B,将B=D-A代入(6)式有：uk+1 = (I-。)1? + D-lb (7)为了保证迭代格式的收敛，需要使迭代矩阵的谱半径p(T) 1.(二)SOR迭代对于(1)式中的线性方程组，将A分裂为：A = D-L-U (8)其中L, U分别为A对应的上三角与下三角矩阵的负值。由Gauss-Siedel迭代格式的中间步骤：D1 声 1 = Luk+1 + Ruk + f (9)引入非零参数3作为松弛因子，则有：Uk+

3、1 = / + 3(小+1 - #)(10)将(9)式代入(10)式，整理之后可以得到SOR迭代格式：uk+1 = (D - a)L)-1(l- 3)。+ a)Ruk + (D -侬工厂1“/ (11) 计算可得最佳松弛因子为：opt2i +Vi-p2(12)(13)(14)其中_p = cosnh(三)Gauss-Siedel 迭代在SOR方法中令松弛因子3 = 1,则有114+1 = #+1则有迭代格式：Duk+1 = Luk+1 + Ruk +/(15)可见Jacobi迭代与Gauss-Siedel迭代的区别在于Jacobi迭代将A分解为D,L+R, Gauss-Siedel 迭代将A分

4、解为D-L, Ro三.程序代码按照上述三种格式分别在matlab中实现，代码如下：(一)Jacobi 迭代说明：u为最终的运算结果，A为系数矩阵，tol为允许的最大迭代步数，steps记录 当前已经经过的迭代步数,eps为规定的计算精度，uO为初值 huiction u,steps = Jacobian iteration( A,uO,b,cps,tol )%u is tlie final result, A is the iteration matrix， tol is tlie pcmiittcd biggest step,%steps shows how many steps of th

5、e iteration,eps is the precision of tlie%resiilt%u0为初值ifnargin5tol=100000;endifnargin=l)fjprintf(invalid iteration:tlie spectral radius is larger tliaii 1);u=uO:k=0;rctiini;end%d=Db;ciT=A*uO-b;%误差的计克采用前后两步迭代之间向量的差，提高计算效率wliile max(abs(eiT)eps&stcps=tol)fjMintitcration exceeds limitjacobian);stepsu=u

6、O;return;endu=iiO;end(二)SOR迭代fiuiction u,steps = SOR_Itcration( A,uO,b,eps,tol )if nargin5 tol=10000;endif nargincps&stcps=tolfpnntfCiteration exceeds limitation,SOR);endend以上代码给出的是一般的SOR迭代方法，具有一般性。实际上针对Poisson方程，mu是可以预先计第 出来的。(三)Gauss-Siedel格式fiuiction u,steps = Gauss_Sicdel Jtcration( A,uO,b,eps,t

7、ol )if nargin5 tol=10000;endif nargincps&stepstol fprintf(fsteps exceed limit,G-S*) return;end end (四)主程序 clc clear %网格剖分 xa=0;xb=l;N1 =64 ;hl =(xb -x a)/Nl;x=xa+l:(Nl-l)*hl;ya=0;yb=l;N2=64;li2=(yb-ya)/N2; y=ya+l:(N2-l)*li2; %gencrate matrix A e2=ones(Nl-l,l);Kl=spdiags(e2,-2*e2,e2,-l,04JSn-lJSri-l)

8、; c3=ones(N2-l,l);Il =spdiags(e3,02-1 JST2-1); A=kion(Kl ,11);A=-A/hlA2;%generate Matrix B c4=oncs(Nl-l,l);I2=spdiags(c4,0Tl-lJSl-l); e5=-2*ones(N2-l ,1);e6=oncs(N2-l,l);K2=spdiags(c6,e5,c6,-l,0,lJSr2-l,N2-l); B=kron(I2JC2);B=B/h2 八 2; %gcnerate g %generate tlie vector of function f f=(x,y)-(2*piA2)

9、*exp(pi*(x+y)*(sin(pi*x)*cos(pi*y)+cos(pi*x)*sin(pi*y); fvec=ones(Nl-l)*(N2-l),l);fori=l:N14for j=l:N2-lf_vec(i-l)*(N2.1)4-j)=f(x(i),y(j);endendX,Y=incshgnd(x,y);uO=zcros0cngtli(fvcc),l);%分别采用三种送代方法求解紧凑格式的近似解ul,stepl=Jacobiaii_Iteration(A+B,uO,fvec);u2,step2=Gauss_Siedel_Itcration(A+BjiO,fvec);u3,st

10、ep3=SOR_Itcration(A+B,uO,fvec);%calciilatc tlic en orlire al= (x ,y) cxp(pi*(x+y)* sin(pi*x).* sin(pi *y);fori=l:Nl-lu_m(i-l)*(N2-l)+l:i*(N2-l)=ii_real(x(i),y);endii_v=ii_inr;ciT_J=inax (ab s(ul -u_v);ciT_G=inax (abs(u2-u_v);en _S=inax (ab s (ii3 -u_v);%月计算结果显示在一张表中res_m=l stcpl err_J;2 stcp2 cn_G;3

11、 stcp3 cit_S;用 rintffMcthod: 1 Jacobiai】2GS3 SOR);fprintethod step error1); resinsol=rcshapc(ul J42-1 mesh(X,Y,sol)四.数值结果针对课本上93页的问题3,程序运行的结果为：表1三种迭代方法性能比较stepserrorJacobian68160. 035788G-S37280. 032068SOR1810. 028627从数值运行结果可以看出，在误差基本相同的情况下，三种迭代方法的收敛速度存在较大差异，其 中Jacobian迭代法需要迭代6846次，约为G-S的二倍，而SOR方法只迭

12、代了 184步就得到了具有 更小误差的解，可见SOR方法在构造上比较复杂，但是实际运行效果却是最佳的。最终迭代结果可 以绘制出图像：图1步长为1/64时利用迭代法的计算结果针对在黑板上布置的问题，其运行结果为:迭代步数误差Jacobian50910. 010989G-S29520. 005419SOR1730. 00047从运算结果上来看，Jacobi迭代方法与GS的步数之间仍然具有:倍关系，SOR迭代的收敛速度更 快并且具有更高的精度。五.计算中出现的问题，解决方法及体会在计算过程中，困扰我的最大问题就是程序运行的效率很低，在最开始一个程序需要运行1到2 分钟的时间，后来在同学的帮助下我们找出了以下问题：矩阵的存储方式：采用稀疏存储spdiags与diag相比会大大提高计算效率，这是因为diag中的0元素参与运算导致程序运行效率低下：误差的计算方式：最初我采用的误差计算方式是用en-Aiib进行计算，这种误差的计克方式导致迟迟无法收敛