数值实验报告Ⅱ
实验名称 | 迭代法求解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。 | ||||||||||||||||||||||||||||||||||||||
三．程序代码 按照上述三种格式分别在matlab中实现，代码如下： （一） Jacobi迭代 说明：u 为最终的运算结果，A 为系数矩阵，tol为允许的最大迭代步数，steps 记录当前已经经过的迭代步数,eps 为规定的计算精度，u0为初值 function [ u,steps ] = Jacobian_Iteration( A,u0,b,eps,tol ) %u is the final result，A is the iteration matrix，tol is the permitted biggest step， %steps shows how many steps of the iteration,eps is the precision of the %result %u0为初值 if nargin<5 tol=100000; end if nargin<4 eps=1e-5; end steps=0; [m,n]=size(A); if m~=n fprintf('the number of columns and rows are unequal'); return end %生成对角矩阵D，采用稀疏存储spdiags而不是diag可以大大节省内存空间并提高计算效率 d=diag(A); e1=ones(n,1); e2=zeros(n,1); D=spdiags([e2 d e2],[-1 0 1],n,n); I=spdiags([e2 e1 e2],[-1 0 1],n,n); %生成迭代矩阵 M=I-(D\A); %{ [vec val]=eig(M); spec_r=max(max(abs(val))); %检验迭代法是否收敛的话需要添加此段代码 if(spec_r>=1) fprintf('invalid iteration:the spectral radius is larger than 1'); u=u0; k=0; return; end %} d=D\b; err=A*u0-b; %误差的计算采用前后两步迭代之间向量的差，提高计算效率 while max(abs(err))>eps&&steps<tol u1=M*u0+d; err=u1-u0; steps=steps+1; u0=u1; end if(steps>=tol) fprintf('iteration exceeds limit,Jacobian'); steps u=u0; return; end u=u0; end （二） SOR迭代 function [ u,steps ] = SOR_Iteration( A,u0,b,eps,tol ) fprintf格式if nargin<5 tol=10000; end if nargin<4 eps=1e-5; end %检查矩阵是否是方阵 [m,n]=size(A); if m~=n fprintf('invalid A:matrix dimentions should agree'); end steps=0; %生成迭代矩阵 d=diag(A); e1=zeros(n,1); e2=ones(n,1); D=spdiags([e1 d e1],[-1 0 1],n,n); I=spdiags([e1 e2 e1],[-1 0 1],n,n); L=-tril(A)+D; R=-triu(A)+D; %这里采取的方法是一般SOR迭代法，需要计算迭代矩阵的特征值，针对特定问题实际上参数是已知的 B=I-D\A; [vec val]=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; %迭代过程 while max(abs(err))>eps&&steps<tol u1=M*u0+d; err=u1-u0; steps=steps+1; u0=u1; end u=u0; if steps>=tol fprintf('iteration exceeds limitation,SOR'); end end 以上代码给出的是一般的SOR迭代方法，具有一般性。实际上针对Poisson方程，mu是可以预先计算出来的。 （三） Gauss-Siedel格式 function [ u,steps ] = Gauss_Siedel_Iteration( A,u0,b,eps,tol ) if nargin<5 tol=10000; end if nargin<4 eps=1e-5; end [m,n]=size(A); steps=0; if(m~=n) fprintf('Invalid A:Matrix dimentions must agree'); u=u0; return; end %generate iteration matrix d=diag(A); e=zeros(n,1); D=spdiags([e d e],[-1 0 1],n,n); L=-tril(A)+D; R=-triu(A)+D; M=(D-L)\R; d=(D-L)\b; err=A*u0-b; while max(abs(err))>eps&&steps<tol u1=M*u0+d; err=u1-u0; steps=steps+1; u0=u1; end u=u0; if steps>tol fprintf('steps exceed limit,G-S') return; end end （四） 主程序 clc clear %网格剖分 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; %generate matrix A e2=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; %generate Matrix B e4=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; %generate g %generate the vector of function f f=@(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); for i=1:N1-1 for j=1:N2-1 f_vec((i-1)*(N2-1)+j)=f(x(i),y(j)); end end [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); %calculate the error u_real=@(x,y)exp(pi*(x+y))*sin(pi*x).*sin(pi*y); for i=1:N1-1 u_m((i-1)*(N2-1)+1:i*(N2-1))=u_real(x(i),y); end u_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=[1 step1 err_J;2 step2 err_G;3 step3 err_S]; fprintf('Method:1-Jacobian,2-GS,3-SOR'); fprintf('Method step error'); res_m sol=reshape(u1,N2-1,N1-1); mesh(X,Y,sol) 四. 数值结果 针对课本上93页的问题3，程序运行的结果为： 表1 三种迭代方法性能比较
从数值运行结果可以看出，在误差基本相同的情况下，三种迭代方法的收敛速度存在较大差异，其中Jacobian迭代法需要迭代6846次，约为G-S的二倍，而SOR方法只迭代了184步就得到了具有更小误差的解，可见SOR方法在构造上比较复杂，但是实际运行效果却是最佳的。最终迭代结果可以绘制出图像： 图1 步长为1/64时利用迭代法的计算结果 针对在黑板上布置的问题，其运行结果为：
从运算结果上来看，Jacobi迭代方法与G-S的步数之间仍然具有二倍关系，SOR迭代的收敛速度更快并且具有更高的精度。 五. 计算中出现的问题，解决方法及体会 在计算过程中，困扰我的最大问题就是程序运行的效率很低，在最开始一个程序需要运行1到2分钟的时间，后来在同学的帮助下我们出了以下问题： 1） 矩阵的存储方式：采用稀疏存储spdiags与diag相比会大大提高计算效率，这是因为diag中的0元素参与运算导致程序运行效率低下； 2） 误差的计算方式：最初我采用的误差计算方式是用err=Au-b进行计算，这种误差的计算方式导致迟迟无法收敛到事先规定好的精度，并且衡量误差是采用err取范数的形式，这又拖慢了计算速度，发现问题之后采取了相邻两步之间迭代向量做差的方式来计算误差，并将误差取为两次迭代之间差最大的分量取绝对值，提升了计算速度； 在发现并解决问题的过程中，不要因为暂时的挫败丧失信心，及时与比自己更优秀的同学交流是更好的学习方式，在交流中汲取知识，共同进步，这也是编程带给我的又一个收获。 | ||||||||||||||||||||||||||||||||||||||
教 师 评 语 | 指导教师： 年 月 日 | |||||||||||||||||||||||||||||||||||||
发表评论