清华大学贾仲孝老师高等数值分析第二次实验.docx
《清华大学贾仲孝老师高等数值分析第二次实验.docx》由会员分享，可在线阅读，更多相关《清华大学贾仲孝老师高等数值分析第二次实验.docx（24页珍藏版）》请在冰豆网上搜索。
清华大学贾仲孝老师高等数值分析第二次实验
高等数值分析第二次实验作业
T1.构造例子特征值全部在右半平面时,观察基本的Arnoldi方法和GMRES方法的数值性态,和相应重新启动算法的收敛性.
Answer:
（1）构造特征值均在右半平面的矩阵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:
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||）'）
得到：
得到两种方法的收敛曲线如上所示，将计算结果整理在下表中：
方法
Arnoldi
GMRES
2.95e-05
3.07e-05
迭代次数
546
526
运行时间（s）
1.564715
92.828966
结果讨论：
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=[];
form=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_restartAMaxiA];
num_restartG=[num_restartGMaxiG];
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||）'）
得到的收敛曲线结果如下图所示：
得到两种方法的收敛曲线如上所示，将计算结果整理在下表中：
方法
ArnoldiM
GMRESM
2.07e-05
2.29e-05
重启次数
26
22
总迭代次数
1300
1100
运行时间（s）
0.844577
0.803231
结果讨论：
1.重启次数随着m的增大而减小，当m增大到一定程度如50之后，减小很缓慢，当m非常大的时候，就过度到了基本算法。
重启的算法如何选择合适的m的因问题而不同。
2.重启算法的总迭代次数（重启次数×m）要多于基本的算法。
这是由于在重启动时，从另一个我们认为更好的初值点（其实不一定好）x0重新开始计算，可能会丢掉一些之前算过的信息。
3.重启算法与基本算法相比，虽然总迭代次数增加了，但是运算速度却能大大提高，运行时间远远小于基本算法（尤其是重启GMRES算法）。
4.一般情况，对于同一个题目，重启的GMRES方法收敛速度要比Arnoldi方法快；
5.总的来说，对题中的矩阵A，四种方法均能稳定地收敛。
T2.对于1中的矩阵,将特征值进行平移,使得实部有正有负,和1的结果进行比较,方法的收敛速度会如何?
基本的Arnoldi算法有无峰点?
若有,基本的GMRES算法相应地会怎样?
Answer:
对A的特征值进行平移，如对S矩阵的实部统一减去一个数s，则小单元矩阵S对应的一对特征值变为
，当矩阵A构造同上题，那么控制s的大小，即控制了实部为负的特征值的个数。
这里将上面的矩阵生成做如下改造：
clearall;
N=500
tol=1e-6;
m=50;
Maxm=50;
A=zeros（2*N）;
s=1.5
fora=1:
N
A（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;
end
U=orth（rand（2*N,2*N））;
A1=U'*A*U;
改变s的值，进而改变实部为负的特征值个数，计算结果整理如下表所示：
实部为负的特征值个数
0
1
5
100
150
200
500
1000
迭代次数
Arnoldi方法
546
638
730
788
462
277
46
16
GMRES方法
526
626
719
757
427
257
45
16
结果讨论：
1.当开始有负半平面的特征值时，个数比较少时（如小于100时），随着个数的增加，基本的Arnoldi和基本的GMRES迭代次数明显变多，收敛速度明显变慢；当负半平面特征值个数继续增加（如大于100时），两者的迭代次数减少，收敛速度明显变快，并且将比第一题的收敛速度快很多，迭代次数少很多。
2.当开始有负半平面的特征值时，个数比较少时（如小于100时），随着个数的增加，Arnoldi出现峰点，峰点个数与其特征值的分布有关系，但整体仍收敛。
这是由于负特征值的存在，使得海森贝格矩阵H发生近似奇异而发生近似中断而引起的。
然而，GMRES的残差始终平稳下降，当Aronldi出现尖峰时，GMRES的残差不变具有最优性。
T3.对1中的例子固定特征值的实部,变化虚部,比较收敛性。
Answer:
首先对T1的代码进行简单修改，实部不变，虚部进行一定的放大（引入系数k）。
每个对角块具有如下形式，对应一对特征向量
上式中的
相应代码如下所示（以Arnoldi为例），取k=0.2,0.5,1,2,5这五种情况。
N=500
A1=zeros（2*N）;
k1=0.2;
fora=1:
N
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;
end
U=orth（rand（2*N,2*N））;
A1=U'*A1*U;
%篇幅所限，此处省去了A2~A5，同理可得。
k5=5;
fora=1:
N
A5（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;
end
U=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'）;
计算结果如下所示：
结果讨论：
1.随着比例因子k的增大，虚部被放大，Arnoldi方法和GMRES方法的收敛速度均减慢，迭代次数增加；
2.同时，随着比例因子k的增大，Arnoldi过程跳动越来越严重，峰点个数增加，发生近似中断的次数增加；而GMRES方法的残差始终保持单调平稳下降。
3.GMRES方法的迭代次数略小于Arnoldi方法的迭代次数，但GMRES的计算时间却远大于Arnoldi方法的计算时间。
T4.当A只有m个不同特征值时,对于大的m和小的m,观察Arnoldi方法和GMRES方法的收敛性.
Answer:
A的生成参考第一次上机作业的第三题的代码，分别构建m=10、50、100、400、800五个矩阵A，代码如下所示：
N=1000
x0=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方法的收敛曲线如下图所示，将相关结果整理在表格中：
不同特征值个数m
10
50
100
400
800
迭代次数
Arnoldi方法
10
30
44
90
97
GMRES方法
10
30
43
88
98
结果讨论：
1.随着不同特征值个数m的增大，Arnoldi方法和GMRES方法的收敛速度均减慢，迭代次数增加。
2.由图可见，Arnoldi方法有小的起伏，而GMRES方法的残差始终单调平稳下降。
同上题一样地，GMRES方法的迭代次数略小于Arnoldi方法的迭代次数。
3.由表格可知，当A只有m个不同的特征值，两种方法至多m步就可以找到精确解；当m较小时，可能需要接近于m步才能找到准确解；而在m较大的时候，算法收敛所需要的迭代次数远小于m。
4.上面的性质与实验1中的Lanczos方法得到的结果类似。
T5.取初始近似解为零向量,右端项b仅由A的m个不同个特征向量的线性组合表示时,Arnoldi方法和GMRES方法的收敛性如何?
Answer:
同本作业第一题构造矩阵A，构造m不同时的右端项b的方法参考第一次上机作业的第四题，如下所示：
N=500
e=1e-6;
m=1000;
A=zeros（2*N）;
fora=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））;
A=U'*A*U;
x0=zeros（2*N,1）;
b=zeros（2*N,1）;
m1=10;
fori=1:
m1
b1=b+rand
（1）*U（:
i）;
end
%篇幅所限，此处省去了b2~b4，同理可得。
m5=800;
fori=1:
m5
b5=b+rand
（1）*U（:
i）;
end
%篇幅所限，此处省去了后续的计算和作图函数，代码同T3一样。
计算得到的基本Arnoldi方法和基本GMRES方法的收敛曲线如下图所示，将相关结果整理在表格中：
不同特征向量个数m
10
50
100
400
800
迭代次数
Arnoldi方法
550
552
548
548
542
GMRES方法
531
529
531
512
524
结果讨论：
1.由结果可知，对于b由m个特征向量组合成的情形，Arnoldi方法和GMRES方法迭代次数和计算时间均无明显变化；
2.由图可见，Arnoldi方法的收敛曲线有小的起伏，而GMRES方法的残差始终单调平稳下降。
同上题一样地，GMRES方法的迭代次数略小于Arnoldi方法的迭代次数，但GMRES的计算时间则大为增加。
3.上面的性质与实验1中的Lanczos方法得到的结果不同了。
附录：
主要算法代码
Arnoldi法Arnoldi.m
function[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:
m
w=A*V（:
j）;
fori=1:
j
H（i,j）=V（:
i）'*w;
w=w-H（i,j）*V（:
i）;
end
H（j+1,j）=norm（w,2）;
ifH（j+1,j）<1e-13
disp（'Arnoldiluckystop'）
x=x0+V（:
1:
j）*ym;
flag=1;
break;
else
V=[V,w/H（j+1,j）];
end
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];
ifrm_norm/norm（b,2）x=x0+V（:
1:
j）*ym;
disp（'ArnoldiConvergeandgetananswer!
'）
flag=1;
break;
end
end
end
GMRES法GMRES.m
function[x,rm,flag]=GMRES（A,b,x0,tol,m）
flag=0;
H=[];
rm=[];
r0=b-A*x0;
nr0=norm（r0,2）;
V=r0/nr0;
x=0;
forj=1:
m
w=A*V（:
j）;
fori=1:
j
H（i,j）=V（:
i）'*w;
w=w-H（i,j）*V（:
i）;
end
H（j+1,j）=norm（w,2）;
ifH（j+1,j）<1e-13
ym=Rm（1:
j,1:
j）\Qm（1:
j,1）*norm（r0,2）;
x=x0+V（:
1:
j）*ym;
flag=1;
break;
else
V=[V,w/H（j+1,j）];
end
Rm=H（1:
j+1,1:
j）;
Qm=eye（j+1,j+1）;
rm_norm1=nr0;
forq=1:
j
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+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;
rm=[rm,rm_norm];
ifrm_norm/norm（b,2）ym=Rm（1:
j,1:
j）\Qm（1:
j,1）*norm（r0,2）;
x=x0+V（:
1:
j）*ym;
flag=1;
break;
end
end
end
m步重启动Arnoldi法ArnoldiM.m
function[x,rm,flag,Maxi]=ArnoldiM（A,b,x0,tol,m,Maxm）
flag=0;
Maxi=1;
rm=[];
whileflag==0
H=[];
r0=b-A*x0;
nr0=norm（r0,2）;
V=r0/nr0;
forj=1:
m
w=A*V（:
j）;
fori=1:
j
H（i,j）=V（:
i）'*w;
w=w-H（i,j）*V（:
i）;
end
H（j+1,j）=norm（w,2）;
ifH（j+1,j）<1e-13
disp（'Arnoldiluckystop'）
x=x0+V（:
1:
j）*ym;
flag=1;
else
V=[V,w/H（j+1,j）];
end
end
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_norm=norm（b-A*x,2）;
rm=[rm,rm_norm];
x=x0+V（:
1:
j）*ym;
x0=x;
ifrm_norm/norm（b,2）disp（'ArnoldiMConvergeandgetananswer!
'）
flag=1;
break;
end
ifMaxi==Maxm
disp（'Arnoldi:
reachmaxsteps!
'）
break;
end
Maxi=Maxi+1;
end
end
m步重启动GMRES法GMRESM.m
function[x,rm,flag,Maxi]=GMRESM（A,b,x0,tol,m,Maxm）
flag=0;
Maxi=1;
rm=[];
whileflag==0
H=[];
r0=b-A*x0;
nr0=norm（r0,2）;
V=r0/nr0;
forj=1:
m
w=A*V（:
j）;
fori=1:
j
H（i,j）=V（:
i）'*w;
w=w-H（i,j）*V（:
i）;
end
H（j+1,j）=norm（w,2）;
ifH（j+1,j）<1e-13
ym=Rm（1:
j,1:
j）\Qm（1:
j,1）*norm（r0,2）;
x=x0+V（:
1:
j）*ym;
flag=1;
break;
else
V=[V,w/H（j+1,j）];
end
Rm=H（1:
j+1,1:
j）;
Qm=eye（j+1,j+1）;
rm_norm1=nr0;
forq=1:
j
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+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
end
ym=Rm（1:
j,1:
j）\Qm（1:
j,1）*norm（r0,2）;
x=x0+V（:
1:
j）*ym;
x0=x;
rm_norm=abs（Qm（j+1,1））*nr0;
rm=[rm,rm_norm];
ifrm_norm/norm（b,2）disp（'GMRESMConvergeandgetananswer!
'）
flag=1;
break;
end
升级会员 