清华大学贾仲孝老师高等数值分析第二次实验.docVIP

WORD完美整理版 范文范例 参考指导 高等数值分析第二次实验作业  T1. 构造例子特征值全部在右半平面时, 观察基本的Arnoldi 方法和GMRES 方法的数值性态, 和相应重新启动算法的收敛性. Answer: 构造特征值均在右半平面的矩阵A： 根据实Schur分解，构造对角矩阵D由n个块形成，每个对角块具有如下形式，对应一对特征值 这样D=diag(S1,S2,S3……Sn)矩阵的特征值均分布在右半平面。生成矩阵A=UTAU，其中U为正交阵，则A矩阵的特征值也均在右半平面。不妨构造A如下所示： 由于选择初值与右端项：x0=zeros(2*N,1);b=ones(2*N,1); 则生成矩阵A的过程代码如下所示： N=500 %生成A为2N阶 A=zeros(2*N); for a=1:N A(2*a-1,2*a-1)=a; A(2*a-1,2*a)=-a; A(2*a,2*a-1)=a; A(2*a,2*a)=a; end U = orth(rand(2*N,2*N)); A1 = U*A*U; 观察基本的Arnoldi和GMRES方法 编写基本的Arnoldi函数与基本GMRES函数，具体代码见附录。 function [x,rm,flag]=Arnoldi(A,b,x0,tol,m) function [x,rm,flag]=GMRES(A,b,x0,tol,m) 输入:A为方程组系数矩阵，b为右端项，x0为初值，tol为停机准则，m为人为限制的最大步数。 输出:x为方程的解，rm为残差向量，flag为解是否收敛的标志。 外程序如下所示： e=1e-6; m=700; tic [xA,rmA,flagA]=Arnoldi(A1,b,x0,e,m); toc tic [xG,rmG,flagG]=GMRES(A1,b,x0,e,m); toc subplot(1,2,1); semilogy(rmA) title(Arnoldiê?á2?ú??) xlabel(μü′ú2?êyk) ylabel(log(||rk||)) subplot(1,2,2); semilogy(rmG) title(GMRESê?á2?ú??) xlabel(μü′ú2?êyk) ylabel(log(||rk||)) 得到： 得到两种方法的收敛曲线如上所示，将计算结果整理在下表中： 方法 Arnoldi GMRES rm 2.95e-05 3.07e-05 迭代次数 546 526 运行时间（s） 1.564715 92.828966 结果讨论： 从图中可以看出，基本的Arnoldi方法经过546步收敛，基本的GMRES方法经过526步收敛，基本的GMRES收敛速度会略快于基本的Arnoldi方法。 从图中可以看出，GMRES方法的的性态较Arnoldi方法更好。Arnoldi方法会有平台和不光滑段，但是由于GMRES具有的残差最优性，GMRES方法平稳地收敛，收敛曲线也更光滑。 观察重新启动的Arnoldi和GMRES方法 在上述两个函数的基础上，编写了重新启动的Arnoldi函数（详情见附录）： function [x,rm,flag,Maxi]=ArnoldiM(A,b,x0,tol,m,Maxm) 输入：m为给定步数,Maxm为人为限制的最大重启次数。 输出：x为方程的解，rm为残差向量，flag为解是否收敛的标志，Maxi为重启次数。 首先编写程序，计算重启次数Maxi与给定步数m的关系，为选取给定步数m给出依据： num_restartA=[]; num_restartG=[]; for m=10:10:150 tic [xG,rmG,flagG,MaxiG]=GMRESM(A,b,x0,tol,m,Maxm); [xA,rmA,flagA,MaxiA]=ArnoldiM(A,b,x0,tol,m,Maxm); num_restartA=[num_restartA MaxiA]; num_restartG=[num_restartG MaxiG]; toc end m1=10:10:150; plot(m1,num_restartA,o-,m1,num_restartG,*-) 从上述结果中看到在m=50之后，重启次数随着给定步长的增加减少很慢，进入饱和。故可以取m=50，最大步数可知在50步以内，运算程序如下所示： m=50; Maxm=50; tic [xA,rmA,flagA,MaxiA]=ArnoldiM(A,b,x0,tol,m,Maxm); toc tic [xG,rmG,flagG,