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

1、xG,rmG,flagG=GMRES（A1,b,x0,e,m）;subplot（1,2,1）;semilogy（rmA）title（Arnoldi）xlabel（kylabel（log（|rk|）subplot（1,2,2）;semilogy（rmG）GMRES得到：得到两种方法的收敛曲线如上所示，将计算结果整理在下表中：方法ArnoldiGMRES2.95e-053.07e-05迭代次数546526运行时间（s）1.56471592.828966结果讨论：1.从图中可以看出，基本的Arnoldi方法经过546步收敛，基本的GMRES方法经过526步收敛，基本的GMRES收敛速度会略快于基本的

2、Arnoldi方法。2.从图中可以看出，GMRES方法的的性态较Arnoldi方法更好。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为重启次数。首先编写程序，计算重启次数Ma

3、xi与给定步数m的关系，为选取给定步数m给出依据：num_restartA=;num_restartG=;for m=10:10:150xG,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;m1=10:150;plot（m1,num_restartA,o-,m1,num_restartG,*- 从上述结果中看到在m=50之后，重启次数随着给

4、定步长的增加减少很慢，进入饱和。故可以取m=50，最大步数可知在50步以内，运算程序如下所示：m=50;Maxm=50;m1=1:MaxiA;m2=1:MaxiG;plot（m1,log10（rmA）,m2,log10（rmG）,得到的收敛曲线结果如下图所示：ArnoldiMGMRESM2.07e-052.29e-05重启次数2622总迭代次数130011000.8445770. 8032311.重启次数随着m的增大而减小，当m增大到一定程度如50之后，减小很缓慢，当m非常大的时候，就过度到了基本算法。重启的算法如何选择合适的m的因问题而不同。2.重启算法的总迭代次数（重启次数m）要多于基本的

5、算法。这是由于在重启动时，从另一个我们认为更好的初值点（其实不一定好）x0重新开始计算，可能会丢掉一些之前算过的信息。3.重启算法与基本算法相比，虽然总迭代次数增加了，但是运算速度却能大大提高，运行时间远远小于基本算法（尤其是重启GMRES算法）。4.一般情况，对于同一个题目，重启的GMRES方法收敛速度要比Arnoldi方法快；5.总的来说，对题中的矩阵A，四种方法均能稳定地收敛。T2. 对于1 中的矩阵, 将特征值进行平移, 使得实部有正有负, 和1 的结果进行比较, 方法的收敛速度会如何? 基本的Arnoldi 算法有无峰点? 若有, 基本的GMRES 算法相应地会怎样?对A的特征值进行

6、平移，如对S矩阵的实部统一减去一个数s，则小单元矩阵S对应的一对特征值变为，当矩阵A构造同上题，那么控制s的大小，即控制了实部为负的特征值的个数。这里将上面的矩阵生成做如下改造：clear all;N=500 tol=1e-6;s=1.5N A（2*a-1,2*a-1）=a-s; A（2*a,2*a）=a-s; 改变s的值，进而改变实部为负的特征值个数，计算结果整理如下表所示：实部为负的特征值个数151002005001000Arnoldi方法6387307884622774616GMRES方法626719757427257451.当开始有负半平面的特征值时，个数比较少时（如小于100时），随

7、着个数的增加，基本的Arnoldi和基本的GMRES迭代次数明显变多，收敛速度明显变慢；当负半平面特征值个数继续增加（如大于100时），两者的迭代次数减少，收敛速度明显变快，并且将比第一题的收敛速度快很多，迭代次数少很多。2.当开始有负半平面的特征值时，个数比较少时（如小于100时），随着个数的增加，Arnoldi出现峰点，峰点个数与其特征值的分布有关系，但整体仍收敛。这是由于负特征值的存在，使得海森贝格矩阵H发生近似奇异而发生近似中断而引起的。然而，GMRES的残差始终平稳下降，当Aronldi出现尖峰时，GMRES的残差不变具有最优性。T3. 对1 中的例子固定特征值的实部, 变化虚部,

8、比较收敛性。首先对T1的代码进行简单修改，实部不变，虚部进行一定的放大（引入系数k）。每个对角块具有如下形式，对应一对特征向量上式中的相应代码如下所示（以Arnoldi为例），取k=0.2,0.5,1,2,5这五种情况。A1=zeros（2*N）;k1=0.2; A1（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;*A1*U;%篇幅所限，此处省去了A2A5，同理可得。k5=5; A5（2*a-1,2*a-1）=a; A5（2*a-1,2*a）=-k5*a; A5（2*a,2*a-1）=k5*a; A

9、5（2*a,2*a）=a;A5 = U*A5*U;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）;size（rmA1,2）;size（rmA2,2）;m3=1:size（rmA3,2）;m4=1:size（rmA4,2）;m5=1:size（rmA5,2）;plo

10、t（m1,log10（rmA1）,m2,log10（rmA2）,m3,log10（rmA3）,m4,log10（rmA4）,m5,log10（rmA5）Arnoldihg1 = legend（k=0.2, k=0.5,k=1k=2k=5）;计算结果如下所示：1.随着比例因子k的增大，虚部被放大，Arnoldi方法和GMRES方法的收敛速度均减慢，迭代次数增加；2.同时，随着比例因子k的增大，Arnoldi过程跳动越来越严重，峰点个数增加，发生近似中断的次数增加；而GMRES方法的残差始终保持单调平稳下降。3.GMRES方法的迭代次数略小于Arnoldi方法的迭代次数，但GMRES的计算时间却远

11、大于Arnoldi方法的计算时间。T4. 当A 只有m 个不同特征值时,对于大的m 和小的m, 观察Arnoldi 方法和GMRES 方法的收敛性.A的生成参考第一次上机作业的第三题的代码，分别构建m=10、50、100、400、800五个矩阵A，代码如下所示：N=1000x0=zeros（N,1）;b=ones（N,1）;m1=10;DIA1=linspace（1,10,m1）;DIA1=DIA1 ones（1,N-m1）;D1=diag（DIA1）;U = orth（rand（N,N）;*D1*U;%篇幅所限，此处省去了A2A4，同理可得。m5=800;DIA5=linspace（1,80

12、0,m4）;DIA5=DIA5 ones（1,N-m4）;D5=diag（DIA5）;*D5*U;%篇幅所限，此处省去了后续的计算和作图函数，代码同T3一样。计算得到的基本Arnoldi方法和基本GMRES方法的收敛曲线如下图所示，将相关结果整理在表格中：不同特征值个数m1050400800304490974388981.随着不同特征值个数m的增大，Arnoldi 方法和 GMRES 方法的收敛速度均减慢，迭代次数增加。2.由图可见，Arnoldi 方法有小的起伏，而 GMRES 方法的残差始终单调平稳下降。同上题一样地，GMRES 方法的迭代次数略小于Arnoldi方法的迭代次数。3.由表格

13、可知，当A只有m个不同的特征值，两种方法至多m步就可以找到精确解；当m较小时，可能需要接近于m步才能找到准确解；而在m较大的时候，算法收敛所需要的迭代次数远小于m。4.上面的性质与实验1中的Lanczos方法得到的结果类似。T5. 取初始近似解为零向量, 右端项b 仅由A 的m 个不同个特征向量的线性组合表示时, Arnoldi 方法和GMRES 方法的收敛性如何?同本作业第一题构造矩阵A，构造m不同时的右端项b的方法参考第一次上机作业的第四题，如下所示：N=500A = Ub=zeros（2*N,1）;for i=1:m1 b1=b+rand（1）*U（:,i）;%篇幅所限，此处省去了b2b

14、4，同理可得。m5 b5=b+rand（1）*U（:不同特征向量个数m5505525485425315295125241.由结果可知，对于b由m个特征向量组合成的情形，Arnoldi方法和GMRES方法迭代次数和计算时间均无明显变化；2.由图可见，Arnoldi 方法的收敛曲线有小的起伏，而 GMRES 方法的残差始终单调平稳下降。同上题一样地，GMRES方法的迭代次数略小于Arnoldi方法的迭代次数，但GMRES的计算时间则大为增加。3.上面的性质与实验1中的Lanczos方法得到的结果不同了。附录：主要算法代码Arnoldi法 Arnoldi.mH=;rm=;V=;x=0;flag=0;

15、r0=b-A*x0;nr0=norm（r0,2）;V=r0/nr0;for j=1:mw=A*V（:,j）;jH（i,j）=V（:,i）*w;w=w-H（i,j）*V（:H（j+1,j）=norm（w,2）;if H（j+1,j）1e-13disp（Arnoldi lucky stopx=x0+V（:,1:j）*ym;flag=1;break;elseV=V,w/H（j+1,j）;gm=nr0*1;zeros（j-1,1）;ym=H（1:j,1:j）gm;rm_norm=H（j+1,j）*abs（zeros（j-1,1）;1*ym）;rm=rm,rm_norm;if rm_norm/norm（

16、b,2）tolArnoldi Converge and get an answer!GMRES法 GMRES.mym=Rm（1:j）Qm（1:j,1）*norm（r0,2）;Rm=H（1:j+1,1:j）;Qm=eye（j+1,j+1）;rm_norm1=nr0;for q=1:ci=Rm（q,q）/（Rm（q,q）2+Rm（q+1,q）2）0.5） ;si=Rm（q+1,q）/（Rm（q,q）2+Rm（q+1,q）2）0.5） ;Rm_temp=Rm;Rm_temp（q ,:）= ci*Rm（q,:）+si*Rm（q+1,:Rm_temp（q+1,:）=-si*Rm（q,:）+ci*Rm（q

17、+1,:Rm=Rm_temp;Qm_temp=Qm;Qm_temp（q ,:）= ci*Qm（q,:）+si*Qm（q+1,:Qm_temp（q+1,:）=-si*Qm（q,:）+ci*Qm（q+1,:Qm=Qm_temp;rm_norm1=abs（si）*rm_norm1;end rm_norm=abs（Qm（j+1,1）*nr0; ym=Rm（1: x=x0+V（:j）*ym ; flag=1; break;m步重启动Arnoldi法 ArnoldiM.mMaxi=1;while flag=0 %rm_norm=norm（b-A*x,2）;x0=x;ArnoldiM Converge and get an answer!if Maxi=Maxm disp（Arnoldi: reach max steps!Maxi=Maxi+1;m步重启动GMRES法 GMRESM.mfunction x,rm,flag,Maxi=GMRESM（A,b,x0,tol,m,Maxm）GMRESM Converge and get an answer!