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

高等数值分析第二次实验作业

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||)')

得到的收敛曲线结果如下图所示：

得到两种方法的收敛曲线如上所示，将计算结果整理在下表中：