下载文档原格式
高等数值分析第二次实验作业
T1.构造例子特征值全部在右半平面时, 观察基本的Arnoldi 方法和GMRES 方法的数值性态, 和相应重新启动算法的收敛性. Answer:
（1） 构造特征值均在右半平面的矩阵A ：
根据实Schur 分解，构造对角矩阵D 由n 个块形成，每个对角块具有如下形式，对应一对特
征值i i i αβ±
i
i i i i S αββα-⎛⎫
= ⎪⎝⎭
这样D=diag(S 1,S 2,S 3……S n )矩阵的特征值均分布在右半平面。生成矩阵A=U T
AU ，其中U 为
正交阵，则A 矩阵的特征值也均在右半平面。不妨构造A 如下所示：
2211112222
/2/2/2/2N N
A n n n n ⨯-⎛⎫
⎪ ⎪ ⎪- ⎪
= ⎪ ⎪ ⎪
- ⎪ ⎪⎝
⎭ 由于选择初值与右端项：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;
（2） 观察基本的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ÊÕÁ²ÇúÏß')
xlabel('µü´ú²½Êýk')
ylabel('log(||rk||)')
subplot(1,2,2);
semilogy(rmG)
title('GMRESÊÕÁ²ÇúÏß')
xlabel('µü´ú²½Êýk')
ylabel('log(||rk||)')
得到：
得到两种方法的收敛曲线如上所示，将计算结果整理在下表中：
结果讨论：
1.从图中可以看出，基本的Arnoldi方法经过546步收敛，基本的GMRES方法经过526
步收敛，基本的GMRES收敛速度会略快于基本的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为重启次数。
首先编写程序，计算重启次数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,MaxiG]=GMRESM(A,b,x0,tol,m,Maxm);
toc
m1=1:MaxiA;
m2=1:MaxiG;
plot(m1,log10(rmA),'o-',m2,log10(rmG),'*-')
title('Á½ÖÖÖØÆôËã·¨µÄÊÕÁ²ÇúÏß')
xlabel('ÖØÆô´ÎÊýk')
ylabel('log(||rk||)')
得到的收敛曲线结果如下图所示：
得到两种方法的收敛曲线如上所示，将计算结果整理在下表中：
结果讨论：
1.重启次数随着m的增大而减小，当m增大到一定程度如50之后，减小很缓慢，当
m非常大的时候，就过度到了基本算法。重启的算法如何选择合适的m的因问题
而不同。
2.重启算法的总迭代次数（重启次数×m）要多于基本的算法。这是由于在重启动时，
从另一个我们认为更好的初值点（其实不一定好）x0重新开始计算，可能会丢掉一
些之前算过的信息。
3.重启算法与基本算法相比，虽然总迭代次数增加了，但是运算速度却能大大提高，
运行时间远远小于基本算法（尤其是重启GMRES算法）。
4.一般情况，对于同一个题目，重启的GMRES方法收敛速度要比Arnoldi方法快；
5.总的来说，对题中的矩阵A，四种方法均能稳定地收敛。
T2. 对于1 中的矩阵, 将特征值进行平移, 使得实部有正有负, 和1 的结果进行比较,
