迭代法实验报告_米瑞琪

1数值实验报告Ⅱ实验名称迭代法求解Poisson方程实验时间2016年4月10日姓名米瑞琪班级信息1303学号1309010304成绩一、实验目的，内容1、理解并掌握三类迭代方法（Jacobi,Gauss-Siedel,SOR）迭代方法的构造原理；2、了解矩阵在matlab中不同的存储方式会导致不同的计算效率，学会采用最高效的方式存储矩阵；3、学会在计算机上实现迭代法，并比较不同方法的效率与误差；二、算法描述（一）Jacobi迭代对于线性方程组：𝐀𝐮=𝐛(𝟏)为了构造迭代格式，可将上式改写为：𝐮=𝐓𝐮+𝐝(𝟐)假定A的对角元均不为0，将A分裂为：𝐀=𝐃−𝐁(𝟑)其中D为：𝐃=𝐝𝐢𝐚𝐠(𝑨𝟏𝟏,𝑨𝟐𝟐,𝑨𝟑𝟑,…,𝑨𝒏𝒏)(𝟒)则(1)式可写为：𝐃𝐮=𝐁𝐮+𝐛(𝟓)形成Jacobi迭代：𝒖𝒌+𝟏=𝑫−𝟏𝑩𝒖𝒌+𝑫−𝟏𝒃(𝟔)为了在迭代格式中不出现B，将B=D-A代入(6)式有：𝒖𝒌+𝟏=(𝐈−𝑫−𝟏𝐀)𝒖𝒌+𝑫−𝟏𝒃(𝟕)为了保证迭代格式的收敛，需要使迭代矩阵的谱半径𝛒(𝐓)𝟏.（二）SOR迭代对于（1）式中的线性方程组，将A分裂为：𝐀=𝐃−𝐋−𝐔(𝟖)其中L，U分别为A对应的上三角与下三角矩阵的负值。由Gauss-Siedel迭代格式的中间步骤：𝑫𝒖𝒌+𝟏̃=𝑳𝒖𝒌+𝟏+𝑹𝒖𝒌+𝒇(𝟗)引入非零参数𝛚作为松弛因子，则有：𝒖𝒌+𝟏=𝒖𝒌+𝝎(𝒖𝒌+𝟏̃−𝒖𝒌)(𝟏𝟎)将（9）式代入（10）式，整理之后可以得到SOR迭代格式：𝒖𝒌+𝟏=(𝑫−𝝎𝑳)−𝟏[(𝟏−𝝎)𝑫+𝝎𝑹]𝒖𝒌+(𝑫−𝝎𝑳)−𝟏𝝎𝒇(𝟏𝟏)计算可得最佳松弛因子为：𝝎𝒐𝒑𝒕=𝟐𝟏+√𝟏−𝝁̅𝟐(𝟏𝟐)其中𝝁̅=𝐜𝐨𝐬𝛑𝐡(𝟏𝟑)（三）Gauss-Siedel迭代在SOR方法中令松弛因子𝛚=𝟏，则有𝐮𝒌+𝟏=𝒖̃𝒌+𝟏(𝟏𝟒)则有迭代格式：𝑫𝒖𝒌+𝟏=𝑳𝒖𝒌+𝟏+𝑹𝒖𝒌+𝒇(𝟏𝟓)可见Jacobi迭代与Gauss-Siedel迭代的区别在于Jacobi迭代将A分解为D，L+R，Gauss-Siedel迭代将A分解为D-L，R。2三．程序代码按照上述三种格式分别在matlab中实现，代码如下：（一）Jacobi迭代说明：u为最终的运算结果，A为系数矩阵，tol为允许的最大迭代步数，steps记录当前已经经过的迭代步数,eps为规定的计算精度，u0为初值function[u,steps]=Jacobian_Iteration(A,u0,b,eps,tol)%uisthefinalresult，Aistheiterationmatrix，tolisthepermittedbiggeststep，%stepsshowshowmanystepsoftheiteration,epsistheprecisionofthe%result%u0为初值ifnargin5tol=100000;endifnargin4eps=1e-5;endsteps=0;[m,n]=size(A);ifm~=nfprintf('thenumberofcolumnsandrowsareunequal');returnend%生成对角矩阵D，采用稀疏存储spdiags而不是diag可以大大节省内存空间并提高计算效率d=diag(A);e1=ones(n,1);e2=zeros(n,1);D=spdiags([e2de2],[-101],n,n);I=spdiags([e2e1e2],[-101],n,n);%生成迭代矩阵M=I-(D\A);%{[vecval]=eig(M);spec_r=max(max(abs(val)));%检验迭代法是否收敛的话需要添加此段代码if(spec_r=1)fprintf('invaliditeration:thespectralradiusislargerthan1');u=u0;k=0;return;end%}d=D\b;err=A*u0-b;%误差的计算采用前后两步迭代之间向量的差，提高计算效率whilemax(abs(err))eps&&stepstolu1=M*u0+d;err=u1-u0;steps=steps+1;u0=u1;endif(steps=tol)fprintf('iterationexceedslimit,Jacobian');stepsu=u0;return;3endu=u0;end（二）SOR迭代function[u,steps]=SOR_Iteration(A,u0,b,eps,tol)ifnargin5tol=10000;endifnargin4eps=1e-5;end%检查矩阵是否是方阵[m,n]=size(A);ifm~=nfprintf('invalidA:matrixdimentionsshouldagree');endsteps=0;%生成迭代矩阵d=diag(A);e1=zeros(n,1);e2=ones(n,1);D=spdiags([e1de1],[-101],n,n);I=spdiags([e1e2e1],[-101],n,n);L=-tril(A)+D;R=-triu(A)+D;%这里采取的方法是一般SOR迭代法，需要计算迭代矩阵的特征值，针对特定问题实际上参数是已知的B=I-D\A;[vecval]=eigs(B);mu=max(max(abs(val)));w=2/(1+(1-mu^2)^0.5);T=D-w*L;M=T\((1-w)*D+w*R);d=T\(w*b);err=A*u0-b;%迭代过程whilemax(abs(err))eps&&stepstolu1=M*u0+d;err=u1-u0;steps=steps+1;u0=u1;endu=u0;ifsteps=tolfprintf('iterationexceedslimitation,SOR');endend以上代码给出的是一般的SOR迭代方法，具有一般性。实际上针对Poisson方程，mu是可以预先计算出来的。（三）Gauss-Siedel格式function[u,steps]=Gauss_Siedel_Iteration(A,u0,b,eps,tol)ifnargin5tol=10000;endifnargin4eps=1e-5;end4[m,n]=size(A);steps=0;if(m~=n)fprintf('InvalidA:Matrixdimentionsmustagree');u=u0;return;end%generateiterationmatrixd=diag(A);e=zeros(n,1);D=spdiags([ede],[-101],n,n);L=-tril(A)+D;R=-triu(A)+D;M=(D-L)\R;d=(D-L)\b;err=A*u0-b;whilemax(abs(err))eps&&stepstolu1=M*u0+d;err=u1-u0;steps=steps+1;u0=u1;endu=u0;ifstepstolfprintf('stepsexceedlimit,G-S')return;endend（四）主程序clcclear%网格剖分xa=0;xb=1;N1=64;h1=(xb-xa)/N1;x=xa+[1:(N1-1)]*h1;ya=0;yb=1;N2=64;h2=(yb-ya)/N2;y=ya+[1:(N2-1)]*h2;%generatematrixAe2=ones(N1-1,1);K1=spdiags([e2,-2*e2,e2],[-1,0,1],N1-1,N1-1);e3=ones(N2-1,1);I1=spdiags(e3,0,N2-1,N2-1);A=kron(K1,I1);A=-A/h1^2;%generateMatrixBe4=ones(N1-1,1);I2=spdiags(e4,0,N1-1,N1-1);e5=-2*ones(N2-1,1);e6=ones(N2-1,1);K2=spdiags([e6,e5,e6],[-1,0,1],N2-1,N2-1);B=kron(I2,K2);B=-B/h2^2;%generateg%generatethevectoroffunctionff=@(x,y)-(2*pi^2)*exp(pi*(x+y))*(sin(pi*x)*cos(pi*y)+cos(pi*x)*sin(pi*y));f_vec=ones((N1-1)*(N2-1),1);5fori=1:N1-1forj=1:N2-1f_vec((i-1)*(N2-1)+j)=f(x(i),y(j));endend[X,Y]=meshgrid(x,y);u0=zeros(length(f_vec),1);%分别采用三种迭代方法求解紧凑格式的近似解[u1,step1]=Jacobian_Iteration(A+B,u0,f_vec);[u2,step2]=Gauss_Siedel_Iteration(A+B,u0,f_vec);[u3,step3]=SOR_Iteration(A+B,u0,f_vec);%calculatetheerroru_real=@(x,y)exp(pi*(x+y))*sin(pi*x).*sin(pi*y);fori=1:N1-1u_m((i-1)*(N2-1)+1:i*(N2-1))=u_real(x(i),y);endu_v=u_m';err_J=max(abs(u1-u_v));err_G=max(abs(u2-u_v));err_S=max(abs(u3-u_v));%将计算结果显示在一张表中res_m=[1step1err_J;2step2err_G;3step3err_S];fprintf('Method:1-Jacobian,2-GS,3-SOR');fprintf('Methodsteperror');res_msol=reshape(u1,N2-1,N1-1);mesh(X,Y,sol)四.数值结果针对课本上93页的问题3，程序运行的结果为：表1三种迭代方法性能比较stepserrorJacobian68460.035788G-S37280.032068SOR1840.028627从数值运行结果可以看出，在误差基本相同的情况下，三种迭代方法的收敛速度存在较大差异，其中Jacobian迭代法需要迭代6846次，约为G-S的二倍，而SOR方法只迭代了184步就得到了具有更小误差的解，可见SOR方法在构造上比较复杂，但是实际运行效果却是最佳的。最终迭代结果可以绘制出图像：图1步长为1/64时利用迭代法的计算结果6针对在黑板上布置的问题，其运行结果为：迭代步数误差Jacobian50910.010989G-S29520.005419SOR1730.00047从运算结果上来看，Jacobi迭代方法与G-S的步数之间仍然具有二倍关系，SOR迭代的收敛速度更快并且具有更高的精度。五.计算中出现的问题，解决方法及体会在计算过程中，困扰我的最大问题就是程序运行的效率很低，在最开始一个程序需要运行1到2分钟的时间，后来在同学的帮助下我们找出了以下问题：1）矩阵的存储方式：采用稀疏存储spdiags与diag相比会大大提高计算效率，这是因为diag中的0元素参与运算导致程序运行效率低下；2）误差的计算方式：最初我采用的误差计算方式是用err=Au-b进行计算，这种误差的计算方式导致迟迟无法收敛到事先规定好的精度，并且衡量误差是采用err取范数的形式，这又拖慢了计算速度，发现问题之后采取了相邻两步之间迭代向量做差的方