《清华大学贾仲孝老师高等数值分析报告第二次实验》由会员分享，可在线阅读，更多相关《清华大学贾仲孝老师高等数值分析报告第二次实验（17页珍藏版）》请在金锄头文库上搜索。
1、高等数值分析第二次实验作业T1构造例子特征值全部在右半平面时，观察基本的Arnoldi方法和GMRES方法的数 值性态, 和相应重新启动算法的收敛性.Answer:(1) 构造特征值均在右半平面的矩阵 A：根据实Schur分解，构造对角矩阵D由n个块形成，每个对角块具有如下形式，对应一对特 征值Q 土 i卩8iI卩iii叮ia丿i这样D二diag(S,S,SS)矩阵的特征值均分布在右半平面。生成矩阵A=UtAU，其中U为1 2 3 n正交阵，则A矩阵的特征值也均在右半平面。不妨构造A如下所示：(1 -11 12 -2A 22n/2 - n 12n / 2 n / 2 丿2 N x2 N由于选择
2、初值与右端项：xO = zeros(2*N,1 );b=ones(2*N,1);则生成矩阵A的过程代码如下所示：N=500 %生成A为2N阶 A=zeros(2*N);for a=1:NA(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;endU = orth(rand(2*N,2*N);A1 = U*A*U;(2) 观察基本的Arnoldi和GMRES方法 编写基本的Arnoldi函数与基本GMRES函数，具体代码见附录。function x,rm,flag=Arnoldi(A,b,x0,tol,m) function
3、x,rm,flag=GMRES(A,b,x0,tol,m)输入:A为方程组系数矩阵，b为右端项，X0为初值，tol为停机准则，m为人为限制的最大 步数。输出：x为方程的解，rm为残差向量，flag为解是否收敛的标志。外程序如下所示e=1e-6;m=700;ticxA,rmA,flagA=Arnoldi(A1,b,x0,e,m); toc ticxG,rmG,flagG=GMRES(A1,b,x0,e,m); tocsubplot(1,2,1); semilogy(rmA) title(ArnoldiE0A2QuIB) xlabel(puzu2Eyk) ylabel(log(|rk|) subp
4、lot(1,2,2);semilogy(rmG) title(GMRESE0A2QuIB) xlabel(pUU2坯Eyk) ylabel(log(|rk|) 得到：迭代歩鳖k迭代步数k得到两种方法的收敛曲线如上所示，将计算结果整理在下表中：方法ArnoldiGMRES2.95e-053.07e-05迭代次数546526运行时间（s）1.56471592.828966结果讨论：1. 从图中可以看出，基本的Arnoldi方法经过546步收敛，基本的GMRES方法经过526 步收敛，基本的GMRES收敛速度会略快于基本的Arnoldi方法。2. 从图中可以看出，GMRES方法的的性态较Arnold
5、i方法更好。Arnoldi方法会有平台 和不光滑段，但是由于GMRES具有的残差最优性，GMRES方法平稳地收敛，收敛 曲线也更光滑。(3) 观察重新启动的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_
6、restartA=;num_restartG=;for m=10:10:150 ticxG,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_resta rtA , o- , m1 , num_resta rtG , *-)从上述结果中看到在m=50之后，重启次数随着
7、给定步长的增加减少很慢，进入饱和。 故可以取m=50，最大步数可知在50步以内，运算程序如下所示：m=50;Maxm=50;ticxA,rmA,flagA,MaxiA=ArnoldiM(A,b,x0,tol,m,Maxm);toc ticxG,rmG,flagG,MaxiG=GMRESM(A,b,x0,tol,m,Maxm); tocm1=1:MaxiA;m2=1:MaxiG;plot(m1,log10(rmA),o-,m2,log10(rmG),*-)title(AOOO0oEapAEOA2QuI)xlabel(O06zIEyk)ylabel(log(|rk|) 得到的收敛曲线结果如下图所示
8、：方法ArnoldiMGMRESM2.07e-052.29e-05重启次数2622总迭代次数13001100运行时间（S）0.8445770.803231结果讨论：1. 重启次数随着m的增大而减小，当m增大到一定程度如50之后，减小很缓慢，当 m 非常大的时候，就过度到了基本算法。重启的算法如何选择合适的 m 的因问题 而不同。2. 重启算法的总迭代次数（重启次数Xm）要多于基本的算法。这是由于在重启动时， 从另一个我们认为更好的初值点（其实不一定好）x0重新开始计算，可能会丢掉一 些之前算过的信息。3. 重启算法与基本算法相比，虽然总迭代次数增加了，但是运算速度却能大大提高， 运行时间远远小
9、于基本算法（尤其是重启GMRES算法）。4. 一般情况，对于同一个题目，重启的GMRES方法收敛速度要比Arnoldi方法快；5. 总的来说，对题中的矩阵A,四种方法均能稳定地收敛。T2. 对于 1 中的矩阵, 将特征值进行平移, 使得实部有正有负, 和 1 的结果进行比较,方法的收敛速度会如何? 基本的 Arnoldi 算法有无峰点? 若有, 基本的 GMRES 算法 相应地会怎样?Answer:对 A 的特征值进行平移，如对 S 矩阵的实部统一减去一个数 s ，则小单元矩阵 S 对应的一 对特征值变为Q -s) + i0，当矩阵A构造同上题，那么控制s的大小，即控制了实部为负.的特征值的个
10、数。- s _厂ii 0.幺-s丿ii这里将上面的矩阵生成做如下改造： clear all;N=500 tol=1e-6;m=50;Maxm=50;A=zeros(2*N);s=1.5 for a=1:NA(2*a-1,2*a-1)=a-s;A(2*a-1,2*a)=-a;A(2*a,2*a-1)=a;A(2*a,2*a)=a-s;endU = orth(rand(2*N,2*N);A1 = U*A*U;改变s的值，进而改变实部为负的特征值个数，计算结果整理如下表所示:z怦心水力生蝕ill嵌弊征灯忆J乜半=山此对吐为1 ；AmHd QMRES7CC实部为负的特征值个数0151001502005
11、001000迭代次数Arnoldi 方法5466387307884622774616GMRES方法5266267197574272574516Qkrb于片半V血的对曲Zj&）A-xdiCWiRES亠一苣I扛諛111独（待址值卷:于左牛面的对数为勺50HM HW I 鶴 ZW ZW W 5WWT这代股家乂100 200 300 400 500 600 700 800 迭代次数 k结果讨论：1. 当开始有负半平面的特征值时，个数比较少时（如小于 100 时），随着个数的增加，基 本的Arnoldi和基本的GMRES迭代次数明显变多，收敛速度明显变慢；当负半平面特征 值个数继续增加（如大于 100
12、时），两者的迭代次数减少，收敛速度明显变快，并且将 比第一题的收敛速度快很多，迭代次数少很多。2. 当开始有负半平面的特征值时,个数比较少时（如小于100时），随着个数的增加,Arnoldi 出现峰点，峰点个数与其特征值的分布有关系，但整体仍收敛。这是由于负特征值的存 在，使得海森贝格矩阵H发生近似奇异而发生近似中断而引起的。然而，GMRES的残差始终平稳下降，当Aronldi出现尖峰时，GMRES的残差不变具有最优性。T3. 对 1 中的例子固定特征值的实部, 变化虚部, 比较收敛性。Answer:首先对T1的代码进行简单修改，实部不变，虚部进行一定的放大（引入系数k）。每个对角块具有如下形
13、式，对应一对特征向量Q + iPiiaiLp.i-Pia.丿i上式中的p = kp.相应代码如下所示(以Arnoldi为例)，取k=0.2,0.5,1,2,5这五种情况。N=500A1=zeros(2*N);k1=0.2;for a=1:NA1(2*a-1,2*a-1)=a;A1(2*a-1,2*a)=-k1*a;A1(2*a,2*a-1)=k1*a;A1(2*a,2*a)=a;endU = orth(rand(2*N,2*N);A1 = U*A1*U;%篇幅所限，此处省夫了 A2A5,同理可得。k5=5;for a=1:NA5(2*a-1,2*a-1)=a;A5(2*a-1,2*a)=-k5
14、*a;A5(2*a,2*a-1)=k5*a;A5(2*a,2*a)=a;endU = orth(rand(2*N,2*N);A5 = U*A5*U;x0=zeros(2*N,1);b=ones(2*N,1);e=1e-6;m=1000;xA1,rmA1,flagA1=Arnoldi(A1,b,x0,e,m);xA2,rmA2,flagA2=Arnoldi(A2,b,x0,e,m);xA3,rmA3,flagA3=Arnoldi(A3,b,x0,e,m);xA4,rmA4,flagA4=Arnoldi(A4,b,x0,e,m);xA5,rmA5,flagA5=Arnoldi(A5,b,x0,e,m); m1=1:size(rmA1,2);m2=1:size(rmA2,2);m3=1:size(rmA3,2);m4=1:size(rmA4,2);log10(rmA4),m5,lm5=1:size(rmA5,2);plot(m1,log10(rmA1),m2,log10(rmA2),m3,log10(rmA3),m4, og10(rmA5)title( uArnoldiEOA2QuI
