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

高等数值分析第二次实验作业高等数值分析第二次实验作业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);fora=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;观察基本的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);toctic[xG,rmG,flagG]=GMRES(A1,b,x0,e,m);tocsubplot(1,2,1);semilogy(rmA)title('ArnoldiÊÕÁ²ÇúÏß')xlabel('µü´ú²½Êýk')ylabel('log(||rk||)')subplot(1,2,2);semilogy(rmG)title('GMRESÊÕÁ²ÇúÏß')xlabel('µü´ú²½Êýk')ylabel('log(||rk||)')得到：得到两种方法的收敛曲线如上所示，将计算结果整理在下中：方法ArnoldiGMRESrm2.95e-053.07e-05迭代次数546526运行时间（s）1.56471592.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=[];form=10:10:150tic[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_restartAMaxiA];num_restartG=[num_restartGMaxiG];tocendm1=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);toctic[xG,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('Á½ÖÖÖØÆôËã·¨µÄÊÕÁ²ÇúÏß')xlabel('ÖØÆô´ÎÊýk')ylabel('log(||rk||)')得到的收敛曲线结果如下图所示：得到两种方法的收敛曲线如上所示，将计算结果整理在下表中：方法ArnoldiMGMRESMrm2.07e-052.29e-05重启次数2622总迭代次数13001100运行时间（s）0.8445770.803231结果讨论：重启次数随着m的增大而减小，当m增大到一定程度如50之后，减小很缓慢，当m非常大的时候，就过度到了基本算法。重启的算法如何选择合适的m的因问题而不同。重启算法的总迭代次数（重启次数×m）要多于基本的算法。这是由于在重启动时，从另一个我们认为更好的初值点（其实不一定好）x0重新开始计算，可能会丢掉一些之前算过的信息。重启算法与基本算法相比，虽然总迭代次数增加了，但是运算速度却能大大提高，运行时间远远小于基本算法（尤其是重启GMRES算法）。一般情况，对于同一个题目，重启的GMRES方法收敛速度要比Arnoldi方法快；总的来说，对题中的矩阵A，四种方法均能稳定地收敛。T2.对于1中的矩阵,将特征值进行平移,使得实部有正有负,和1的结果进行比较,方法的收敛速度会如何?基本的Arnoldi算法有无峰点?若有,基本的GMRES算法相应地会怎样?Answer:对A的特征值进行平移，如对S矩阵的实部统一减去一个数s，则小单元矩阵S对应的一对特征值变为，当矩阵A构造同上题，那么控制s的大小，即控制了实部为负的特征值的个数。这里将上面的矩阵生成做如下改造：clearall;N=500tol=1e-6;m=50;Maxm=50;A=zeros(2*N);s=1.5fora=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的值，进而改变实部为负的特征值个数，计算结果整理如下表所示：0100200300400500600700800迭代次数k-5-4-3-2-10123log(||rk||)收敛曲线(特征值位于左半平面的对数为100)ArnoldiGMRES实部为负的特征值个数0151001502005001000迭代次数Arnoldi方法5466387307884622774616GMRES方法5266267197574272574516结果讨论：当开始有负半平面的特征值时，个数比较少时（如小于100时），随着个数的增加，基本的Arnoldi和基本的GMRES迭代次数明显变多，收敛速度明显变慢；当负半平面特征值个数继续增加（如大于100时），两者的迭代次数减少，收敛速度明显变快，并且将比第一题的收敛速度快很多，迭代次数少很多。当开始有负半平面的特征值时，个数比较少时（如小于100时），随着个数的增加，Arnoldi出现峰点，峰点个数与其特征值的分布有关系，但整体仍收敛。这是由于负特征值的存在，使得海森贝格矩阵H发生近似奇异而发生近似中断而引起的。然而，GMRES的残差始终平稳下降，当Aronldi出现尖峰时，GMRES的残差不变具有最优性。T3.对1中的例子固定特征值的实部,变化虚部,比较收敛性。Answer:首先对T1的代码进行简单修改，实部不变，虚部进行一定的放大（引入系数k）。每个对角块具有如下形式，对应一对特征向量上式中的相应代码如下所示（以Arnoldi为例），取k=0.2,0.5,1,2,5这五种情况。N=500A1=zeros(2*N);k1=0.2;fora=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;%篇幅所限，此处省去了A2~A5，同理可得。k5=5;fora=1:NA5(2*a-1,2*a-1)=a;A5(2*a-1,2*a)=-k5*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);m5=1:size(rmA5,2);plot(m1,log10(rmA1),m2,log10(rmA2),m3,log10(rmA3),m4,log10(rmA4),m5,log10(rmA5))title('»ù±¾Arnoldi·½·¨ÊÕÁ²ÇúÏß')xlabel('µü´ú´ÎÊýk')ylabel('log(||rk||)')hg1=legend('k=0.2','k=0.5','k=1','k=2','k=5');计算结果如下所示：结果讨论：随着比例因子k的增大，虚部被放大，Arnoldi方法和GMRES方法的收敛速度均减慢，迭代次数增加；同时，随着比例因子k的增大，Arnoldi过程跳动越来越严重，峰点个数增加，发生近似中断的次数增加；而GMRES方法的残差始终保持单调平稳下降。GMRES方法的迭代次数略小于Arnoldi方法的迭代次数，但GMRES的计算时间却远大于Arnoldi方法的计算时间。T4.当A只有m个不同特征值时,对于大的m和小的m,观察Arnoldi方法和GMRES方法的收敛性.Answer:A的生成参考第一次上机作业的第三题的代码，分别构建m=10、50、100、400、800五个矩阵A，代码如下所示：N=1000x0=zeros(N,1);b=ones(N,1);e=1e-6;m1=10;DIA1=linspace(1,10,m1);DIA1=[DIA1ones(1,N-m1)];D1=diag(DIA1);U=orth(rand(N,N));A1=U'*D1*U;%篇幅所限，此处省去了A2~A4，同理可得。m5=800;DIA5=linspace(1,800,m4);DIA5=[DIA5ones(1,N-m4)];D5=diag(DIA5);A5=U'*D5*U;%篇幅所限，此处省去了后续的计算和作图函数，代码同T3一样。计算得到的基本Arnoldi方法和基本GMRES方法的收敛曲线如下图所示，将相关结果整理在表格中：不同特征值个数m1050100400800迭代次数Arnoldi方法1030449097GMRES方法1030438898结果讨论：随着不同特征值个数m的增大，Arnoldi方法和GMRES方法的收敛速度均减慢，迭代次数增加。由图可见，Arnoldi方法有小的起伏，而GMRES方法的残差始终单调平稳下降。同上题一样地，GMRES方法的迭代次数略小于Arnoldi方法的迭代次数。由表格可知，当A只有m个不同的特征值，两种方法至多m步就可以找到精确解；当m较小时，可能需要接近于m步才能找到准确解；而在m较大的时候，算法收敛所需要的迭代次数远小于m。上面的性质与实验1中的Lanczos方法得到的结果类似。T5.取初始近似解为零向量,右端项b仅由A的m个不同个特征向量的线性组合表示时,Arnoldi方法和GMRES方法的收敛性如何?Answer:同本作业第一题构造矩阵A，构造m不同时的右端项b的方法参考第一次上机作业的第四题，如下所示：N=500e=1e-6;m=1000;A=zeros(2*N);fora=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));A=U'*A*U;x0=zeros(2*N,1);b=zeros(2*N,1);m1=10;fori=1:m1b1=b+rand(1)*U(:,i);end%篇幅所限，此处省去了b2~b4，同理可得。m5=800;fori=1:m5b5=b+rand(1)*U(:,i);end%篇幅所限，此处省去了后续的计算和作图函数，代码同T3一样。计算得到的基本Arnoldi方法和基本GMRES方法的收敛曲线如下图所示，将相关结果整理在表格中：不同特征向量个数m1050100400800迭代次数Arnoldi方法550552548548542GMRES方法531529531512524结果讨论：由结果可知，对于b由m个特征向量组合成的情形，Arnoldi方法和GMRES方法迭代次数和计算时间均无明显变化；由图可见，Arnoldi方法的收敛曲线有小的起伏，而GMRES方法的残差始终单调平稳下降。同上题一样地，GMRES方法的迭代次数略小于Arnoldi方法的迭代次数，但GMRES的计算时间则大为增加。上面的性质与实验1中的Lanczos方法得到的结果不同了。附录：主要算法代码Arnoldi法Arnoldi.mfunction[x,rm,flag]=Arnoldi(A,b,x0,tol,m)H=[];rm=[];V=[];x=0;flag=0;r0=b-A*x0;nr0=norm(r0,2);V=r0/nr0;forj=1:mw=A*V(:,j);fori=1:jH(i,j)=V(:,i)'*w;w=w-H(i,j)*V(:,i);endH(j+1,j)=norm(w,2);ifH(j+1,j)<1e-13disp('Arnoldiluckystop')x=x0+V(:,1:j)*ym;flag=1;break;elseV=[V,w/H(j+1,j)];endgm=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];ifrm_norm/norm(b,2)