清华大学贾仲孝老师高等数值分析第二次实验Word格式.docx
《清华大学贾仲孝老师高等数值分析第二次实验Word格式.docx》由会员分享,可在线阅读,更多相关《清华大学贾仲孝老师高等数值分析第二次实验Word格式.docx(24页珍藏版)》请在冰豆网上搜索。
[xG,rmG,flagG]=GMRES(A1,b,x0,e,m);
subplot(1,2,1);
semilogy(rmA)
title('
ArnoldiÊ
Õ
Á
²
Ç
ú
Ï
ß
'
)
xlabel('
µ
ü
´
½
Ê
ý
k'
ylabel('
log(||rk||)'
subplot(1,2,2);
semilogy(rmG)
GMRESÊ
得到:
得到两种方法的收敛曲线如上所示,将计算结果整理在下表中:
方法
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
[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];
m1=10:
150;
plot(m1,num_restartA,'
o-'
m1,num_restartG,'
*-'
从上述结果中看到在m=50之后,重启次数随着给定步长的增加减少很慢,进入饱和。
故可以取m=50,最大步数可知在50步以内,运算程序如下所示:
m=50;
Maxm=50;
m1=1:
MaxiA;
m2=1:
MaxiG;
plot(m1,log10(rmA),'
m2,log10(rmG),'
Ö
Ø
Æ
ô
Ë
ã
·
¨
Ä
Î
得到的收敛曲线结果如下图所示:
ArnoldiM
GMRESM
2.07e-05
2.29e-05
重启次数
26
22
总迭代次数
1300
1100
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算法相应地会怎样?
对A的特征值进行平移,如对S矩阵的实部统一减去一个数s,则小单元矩阵S对应的一对特征值变为
,当矩阵A构造同上题,那么控制s的大小,即控制了实部为负的特征值的个数。
这里将上面的矩阵生成做如下改造:
clearall;
N=500
tol=1e-6;
s=1.5
N
A(2*a-1,2*a-1)=a-s;
A(2*a,2*a)=a-s;
改变s的值,进而改变实部为负的特征值个数,计算结果整理如下表所示:
实部为负的特征值个数
1
5
100
200
500
1000
Arnoldi方法
638
730
788
462
277
46
16
GMRES方法
626
719
757
427
257
45
1.当开始有负半平面的特征值时,个数比较少时(如小于100时),随着个数的增加,基本的Arnoldi和基本的GMRES迭代次数明显变多,收敛速度明显变慢;
当负半平面特征值个数继续增加(如大于100时),两者的迭代次数减少,收敛速度明显变快,并且将比第一题的收敛速度快很多,迭代次数少很多。
2.当开始有负半平面的特征值时,个数比较少时(如小于100时),随着个数的增加,Arnoldi出现峰点,峰点个数与其特征值的分布有关系,但整体仍收敛。
这是由于负特征值的存在,使得海森贝格矩阵H发生近似奇异而发生近似中断而引起的。
然而,GMRES的残差始终平稳下降,当Aronldi出现尖峰时,GMRES的残差不变具有最优性。
T3.对1中的例子固定特征值的实部,变化虚部,比较收敛性。
首先对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;
%篇幅所限,此处省去了A2~A5,同理可得。
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;
A5(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);
plot(m1,log10(rmA1),m2,log10(rmA2),m3,log10(rmA3),m4,log10(rmA4),m5,log10(rmA5))
»
ù
±
¾
Arnoldi·
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方法的收敛性.
A的生成参考第一次上机作业的第三题的代码,分别构建m=10、50、100、400、800五个矩阵A,代码如下所示:
N=1000
x0=zeros(N,1);
b=ones(N,1);
m1=10;
DIA1=linspace(1,10,m1);
DIA1=[DIA1ones(1,N-m1)];
D1=diag(DIA1);
U=orth(rand(N,N));
*D1*U;
%篇幅所限,此处省去了A2~A4,同理可得。
m5=800;
DIA5=linspace(1,800,m4);
DIA5=[DIA5ones(1,N-m4)];
D5=diag(DIA5);
*D5*U;
%篇幅所限,此处省去了后续的计算和作图函数,代码同T3一样。
计算得到的基本Arnoldi方法和基本GMRES方法的收敛曲线如下图所示,将相关结果整理在表格中:
不同特征值个数m
10
50
400
800
30
44
90
97
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方法的收敛性如何?
同本作业第一题构造矩阵A,构造m不同时的右端项b的方法参考第一次上机作业的第四题,如下所示:
N=500
A=U'
b=zeros(2*N,1);
fori=1:
m1
b1=b+rand
(1)*U(:
i);
%篇幅所限,此处省去了b2~b4,同理可得。
m5
b5=b+rand
(1)*U(:
不同特征向量个数m
550
552
548
542
531
529
512
524
1.由结果可知,对于b由m个特征向量组合成的情形,Arnoldi方法和GMRES方法迭代次数和计算时间均无明显变化;
2.由图可见,Arnoldi方法的收敛曲线有小的起伏,而GMRES方法的残差始终单调平稳下降。
同上题一样地,GMRES方法的迭代次数略小于Arnoldi方法的迭代次数,但GMRES的计算时间则大为增加。
3.上面的性质与实验1中的Lanczos方法得到的结果不同了。
附录:
主要算法代码
Arnoldi法Arnoldi.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);
j
H(i,j)=V(:
i)'
*w;
w=w-H(i,j)*V(:
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)];
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)<
tol
ArnoldiConvergeandgetananswer!
GMRES法GMRES.m
ym=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;
forq=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+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.m
Maxi=1;
whileflag==0
%rm_norm=norm(b-A*x,2);
x0=x;
ArnoldiMConvergeandgetananswer!
ifMaxi==Maxm
disp('
Arnoldi:
reachmaxsteps!
Maxi=Maxi+1;
m步重启动GMRES法GMRESM.m
function[x,rm,flag,Maxi]=GMRESM(A,b,x0,tol,m,Maxm)
GMRESMConvergeandgetananswer!