高等数值分析作业第二次实验.docx
《高等数值分析作业第二次实验.docx》由会员分享，可在线阅读，更多相关《高等数值分析作业第二次实验.docx（24页珍藏版）》请在冰豆网上搜索。
高等数值分析作业第二次实验
高等数值分析第二次实验作业
注：
矩阵阶数均为1000阶
T1.构造例子特征值全部在右半平面时，观察基本的Arnoldi方法和GMRES方法的数值性态，和相应重新启动算法的收敛性。
Answer：
Ø关于特征值均在右半平面的矩阵A：
首先构造对角矩阵D，其中D矩阵的对角是由以下形式的2×2的子矩阵组成其对角元：
此时，S矩阵的特征值分别为a+bi和a-bi，这样D=diag（S1,S2,S3……Sn）矩阵的特征值均分布在右半平面。
另矩阵A=QTAQ，则A矩阵的特征值也均在右半平面，生成矩阵的过程代码如下所示：
N=500%生成的矩阵为2N阶
A=zeros（2*N）;
%delta控制特征值分布的密集程度，即控制条件数大小，delta越大条件数越大
delta=0.1;
%构造D矩阵
forj=1:
N
A（2*j-1,2*j-1）=N+j*delta;
A（2*j-1,2*j）=-N-j*delta;
A（2*j,2*j-1）=N+j*delta;
A（2*j,2*j）=N+j*delta;
end
U=orth（rand（2*N,2*N））;
A=U'*A*U;
Ø首先进行观察基本的Arnoldi方法和GMRES方法的数值性态，考虑N=500，即矩阵为1000阶，取X0=zeros（N,1），b=ones（N,1），收敛准则为e=10-6;
Arnoldi方法函数如下：
function[xm,error,num]=Arnoldi（A,x0,b,e）
N=size（A,1）;
r=b-A*x0;
belta=norm（r）;
v=r/belta;
V=v;
j=0;
whilenorm（r）>e&jj=j+1;
num（j）=j;
w=A*V（1:
N,j）;
fori=1:
j
h（i,j）=V（1:
N,i）'*w;
w=w-h（i,j）*V（1:
N,i）;
end
h（j+1,j）=norm（w）;
ifh（j+1,j）~=0
v=w/h（j+1,j）;
V=[Vv];
end
e1=zeros（j,1）;
e1
（1）=1;
be1=belta*e1;
try
[L,U,S]=lu（h（1:
j,1:
j））;
be1=S*be1;
lym=LX（L,be1）;
ym=UX（U,lym）;
xm=x0+V（1:
N,1:
j）*ym;
r=b-A*xm;
end
error（j）=log10（norm（r））;
end
end
GMRES方法的函数如下：
function[xm,error,num]=GMRES（A,x0,b,e）
%ARNOLDISummaryofthisfunctiongoeshere
%Detailedexplanationgoeshere
N=size（A,1）;
r=b-A*x0;
belta=norm（r）;
v=r/belta;
V=v;
j=0;
er=1000;
whileer>e&jj=j+1;
num（j）=j;
w=A*V（1:
N,j）;
fori=1:
j
h（i,j）=V（1:
N,i）'*w;
w=w-h（i,j）*V（1:
N,i）;
end
h（j+1,j）=norm（w）;
ifh（j+1,j）~=0
v=w/h（j+1,j）;
V=[Vv];
end
try
[Q,R]=qr（h（1:
j+1,1:
j））;
gk=Q'*belta*eye（j+1,1）;
ym=minresYK（R,gk）;
xm=x0+V（1:
N,1:
j）*ym;
er=gk（j+1）;
end
error（j）=log10（er）;
end
end
其中UX，LX和minresYK是自己编写的求解（上、下）三角矩阵的函数，代码如下：
functionyk=LX（TK,b）
n=size（TK,1）;
yk=zeros（n,1）;
yk
（1）=b
（1）/TK（1,1）;
fori=2:
n
yk（i）=b（i）;
forj=1:
i-1
yk（i）=yk（i）-TK（i,j）*yk（j）;
end
yk（i）=yk（i）/TK（i,i）;
end
end
functionyk=UX（TK,b）
n=size（TK,2）;
yk=zeros（n,1）;
yk（n）=b（n）/TK（n,n）;
fori=1:
n-1
k=n-i;
yk（k）=b（k）;
forj=k+1:
n
yk（k）=yk（k）-TK（k,j）*yk（j）;
end
yk（k）=yk（k）/TK（k,k）;
end
end
functionyk=minresYK（TK,b）
n=size（TK,2）;
yk=zeros（n,1）;
yk（n）=b（n）/TK（n,n）;
fori=1:
n-1
k=n-i;
yk（k）=b（k）;
forj=k+1:
n
yk（k）=yk（k）-TK（k,j）*yk（j）;
end
yk（k）=yk（k）/TK（k,k）;
end
end
两种方法的计算结果如下所示：
Delta=1即条件数较大
Delta=0.1即条件数较小
可以看到delta=0.1时两种方法收敛都很快，delta=1时，两种方法略有差别，GMRES方法收敛速度略快一些，下面就去delta=1的情况进行第一题的讨论。
delta=1时，计算结果如下：
方法
Arnoldi
GMRES方法
-6.4922
-6.0110
迭代次数
25
26
运行时间（s）
0.109216
0.135968
可以看到两种方法的迭代次数差不多，有曲线可以看到GMRES方法的的性态较Arnoldi方法良好，Arnoldi方法会有平台段，但是GMRES方法的则平稳地收敛，这也是最后迭代次数GMRES较Arnoldi方法少一次的原因。
Ø对于重启的Arnoldi方法和GMRES方法代码如下所示：
A、首先对之前的标准Arnoldi方法和GMRES方法进行修改，代码如下（m为间隔m步重启）：
function[xm,error]=Arnoldi_m（A,x0,b,e,m）
N=size（A,1）;
r=b-A*x0;
belta=norm（r）;
v=r/belta;
V=v;
j=0;
whileabs（norm（r））>e&jj=j+1;
w=A*V（1:
N,j）;
fori=1:
j
h（i,j）=V（1:
N,i）'*w;
w=w-h（i,j）*V（1:
N,i）;
end
h（j+1,j）=norm（w）;
ifh（j+1,j）~=0
v=w/h（j+1,j）;
V=[Vv];
end
e1=zeros（j,1）;
e1
（1）=1;
be1=belta*e1;
try
[L,U,S]=lu（h（1:
j,1:
j））;
be1=S*be1;
lym=LX（L,be1）;
ym=UX（U,lym）;
xm=x0+V（1:
N,1:
j）*ym;
r=b-A*xm;
end
end
error=abs（norm（r））;
end
function[xm,error]=GMRES_m（A,x0,b,e,m）
N=size（A,1）;
r=b-A*x0;
belta=norm（r）;
v=r/belta;
V=v;
j=0;
er=1000;
whileabs（er）>e&jj=j+1;
w=A*V（1:
N,j）;
fori=1:
j
h（i,j）=V（1:
N,i）'*w;
w=w-h（i,j）*V（1:
N,i）;
end
h（j+1,j）=norm（w）;
ifh（j+1,j）~=0
v=w/h（j+1,j）;
V=[Vv];
end
try
[Q,R]=qr（h（1:
j+1,1:
j））;
gk=Q'*belta*eye（j+1,1）;
ym=minresYK（R,gk）;
xm=x0+V（1:
N,1:
j）*ym;
er=gk（j+1）;
end
error=abs（er）;
end
end
B、调用上面更改后的代码，编写重启的算法代码如下：
function[xm,err,num]=Arnoldim（A,x0,b,e,m）
N=size（A,1）;
jmax=N/m;
j=0;
xm=x0;
error=1000;
whileje
j=j+1;
[xm,error]=Arnoldi_m（A,xm,b,e,m）;
num（j）=j;
err（j）=log10（error）;
end
j
log10（error）
end
function[xm,err,num]=GMRESm（A,x0,b,e,m）
N=size（A,1）;
jmax=N/m;
j=0;
xm=x0;
error=1000;
whileje
j=j+1;
[xm,error]=GMRES_m（A,xm,b,e,m）;
num（j）=j;
err（j）=log10（error）;
end
j
log10（error）
end
计算结果如下所示（计算时取m=5，即每隔5步进行算法的重启）：
delta=1m=5（每隔5步重启）
delta=1m=10（每隔10步重启）
delta=10m=5（每隔5步重启）
delta=10m=10（每隔10步重启）
当delta=100m=5（每隔5步重启），结果如下：
从结果可以看出：
1.重启算法的总迭代次数（重启次数×m）要多于基本的算法；
2.重启的Arnoldi方法收敛速度要比GMRES方法快；
3.对于这个问题，重启算法也可以稳定地收敛，及每次重启后的x0都能有所改善；
4.当条件数增大时，明显重启次数会增加；
5.总的来说，对题中的矩阵A，四种方法均能稳定地收敛。
T2.对于1中的矩阵，将特征值进行平移，使得实部有正有负，和1的结果进行比较，方法的收敛速度会如何？
基本的Arnoldi算法有无峰点？
若有，基本的GMRES算法相应地会怎样？
Answer：
Ø对A的特征值进行平移，由于只要变换A矩阵的实部，这里将上面的矩阵生成做如下改造：
A=zeros（2*N）;
delta=0.1;
s=1.5%s用于控制实部为负的特征值的个数，s越大，实部为负的特征值个数越多。
forj=1:
N
A（2*j-1,2*j-1）=-s*delta+j*delta;
A（2*j-1,2*j）=s*delta-j*delta;
A（2*j,2*j-1）=-s*delta+j*delta;
A（2*j,2*j）=-s*delta+j*delta;
end
U=orth（rand（2*N,2*N））;
A=U'*A*U;
Ø基本算法的结果如下所示：
2个实部为负的特征值
10个实部为负的特征值
100个实部为负的特征值
200个实部为负的特征值
500个实部为负的特征值
900个实部为负的特征值
990个实部为负的特征值
998个实部为负的特征值
1000个实部为负的特征值
实部为负的特征值个数
2
10
100
200
500
Arnoldi方法
-6.0416
-6.0025
-11.8957
-12.2160
-12.2839
迭代次数
734
883
1000
1000
1000
运行时间（s）
17.517751
31.238844
41.911802
41.605778
41.496951
GMRES方法
-6.0008
-6.0702
-12.4162
-12.6358
-12.9732
迭代次数
726
880
1000
1000
1000
运行时间（s）
15.142969
27.069932
37.133932
36.937593
36.433673
实部为负的特征值个数
900
990
998
1000
Arnoldi方法
-12.0507
-6.0954
-6.0192
-6.0251
迭代次数
1000
882
738
656
运行时间（s）
42.070828
28.833792
17.203467
12.627750
GMRES方法
-12.5621
-6.0702
-6.0105
-6.0043
迭代次数
1000
878
730
646
运行时间（s）
36.308397
24.999022
16.011876
10.850091
Ø重启算法的计算结果如下所示：
2个实部为负的特征值，每10步重启
2个实部为负的特征值，每100步重启
10个实部为负的特征值，每10步重启
10个实部为负的特征值，每100步重启
100个实部为负的特征值，每10步重启
100个实部为负的特征值，每100步重启
500个实部为负的特征值，每10步重启
500个实部为负的特征值，每100步重启
998个实部为负的特征值，每10步重启
500个实部为负的特征值，每100步重启
Ø结论：
1.对于上面的矩阵，基本的Arnoldi方法与GMRES方法均能收敛，通知GMRES方法收敛速度相对较快，而重启的算法则均不收敛；
2.可以看到单特征值均匀地分布在左右两边平面时（即有一半的特征值实部为负，另一半实部为正），问题的病态程度最高；
3.问题越病态，重启的Arnoldi算法得到的结果越不好，而重启的GMRES对于结果基本就没有改进，这两种方法均不收敛；
4.基本的Arnoldi算法有峰点，这是由于实部为负的特征值的存在，使得H矩阵发生近似奇异导致发生近似中断而引起的，而GMRES算法则相对稳定，当然这两种算法都能收敛。
5.Arnoldi方法的峰点个数与其特征值的分布有关系，特征在左右半平面分布越均匀，峰点个数越多，当实部为负数的特征值小于N/2（N为矩阵阶数），峰点个数等于这类特征值个数，当实部为负数的特征值大于N/2（N为矩阵阶数），峰点个数等于N减这类特征值个数。
6.两种基本算法收敛速度在前期均很慢，但是到迭代后期收敛速度突增，同时可以看到，当问题相当病态时（例如上面100个、200个、500个实部为负的特征值对应的情况），两种方法均是最后一步才一下子收敛。
T3.对1中的例子固定特征值的实部，变化虚部，比较收敛性。
Answer：
Ø首先对T1的代码进行简单修改，实部不变，虚部进行一定的放大（引入系数α），相应代码如下：
N=500
A=zeros（2*N）;
x0=zeros（2*N,1）;
b=ones（2*N,1）;
e=10^（-6）;
delta=1;
alpha=1;
forj=1:
N
A（2*j-1,2*j-1）=N+j*delta;
A（2*j-1,2*j）=-N-alpha*j*delta;
A（2*j,2*j-1）=N+alpha*j*delta;
A（2*j,2*j）=N+j*delta;
end
U=orth（rand（2*N,2*N））;
A=U'*A*U;
计算结果如下所示，取delta=1，α分别为1、3、5，取X0=zeros（N,1），b=ones（N,1），收敛准则为e=10-6：
Arnoldi方法
GMRES
重启Arnoldi方法，每5步重启
重启GMRES方法，每5步重启
重启Arnoldi方法，每10步重启
重启GMRES方法，每10步重启
重启GMRES方法，每10步重启：
Ø结论：
1.随着虚部被放大，四种算法的收敛速度均越来越慢；
2.虽然速度变慢但是，两种基本的算法均能收敛，
与迭代次数基本成线性关系，而重启GMRES算法也都能收敛，可以看到即使是α为100的情况，重启GMRES算法每一次重启得到的
随着重启次数基本线性减小，即重启后的x0会有所改善；
3.对于重启Arnoldi算法，他是不稳定的，当α过大，而重启间隔的步数较小时，可能发生不收敛的情况，例如上面的α=5，m-5，重启后的结果可能越来越差，并没有得到改善；
T4.当A只有m个不同特征值时，对于大的m和小的m，观察Arnoldi方法和GMRES方法的收敛性.
Answer：
A的生成参考第一次实验作业的代码，如下所示：
N=1000
m=10;
DIA=linspace（1,1000,m）;
VEC=zeros（N,1）;
k=N/m;
i=1;
ii=0;
whilei<=m
forj=1:
k
VEC（ii+j）=DIA（i）;
end
ii=ii+k;
i=i+1;
end
D=diag（VEC）;
U=orth（rand（N,N））;
A=U'*D*U;
对构造好的A进行计算，在计算时，取X0=zeros（N,1）;b=ones（N,1）;e=10-6，计算结果如下：
Arnoldi方法
GMRES方法
相同的特征值个数
10
20
50
100
500
1000
Arnoldi方法
-10.6206
-10.9011
-6.0109
-6.2994
-6.0594
-6.0024
迭代次数
10
20
42
62
128
165
运行时间（s）
0.021357
0.038060
0.082905
0.132112
0.370289
0.584533
GMRES方法
-10.6400
-11.0604
-6.0084
-6.0891
-6.0663
-6.0048
迭代次数
10
20
42
61
125
162
运行时间（s）
0.013224
0.023461
0.055279
0.093193
0.270723
0.447675
显然，上面构造出来的矩阵的条件数是相同的，得到以下几个结论：
1.对于基本的Arnoldi方法和GMRES方法，计算结果总是收敛的，GMRES方法收敛性较Arnoldi方法好，当m较大时，GMRES总能比Arnoldi方法提前几步收敛；
2.如果A只有m个不同的特征值，两种方法至多m步就可以找到精确解；
3.在m较大的时候，算法收敛所需要的迭代次数远小于m，当m较小时，可能需要接近于m步才能找到准确解；
4.由上面计算可以看到在m值较小时，Arnoldi方法计算出来的第k步的误差
会有一个在计算时会有先增大后减小的过程，最后收敛，而GMRES方法则没有这一性质，结果单调地收敛到精确解；
5.上面的性质与实验1中的Lanczos方法得到的结果类似。
T5.取初始值近似解为零向量，右端项b仅由A的m个不同特征向量的线性组合表示时，Arnoldi方法和GMRES方法的收敛性如何？
Answer：
Ø同上面那题一样，参考第一次实验的矩阵A生成的方法，我们可以知道U的每一个行向量均为对应A的特征值（D的某一个元素）的特征向量，因此构造m不同时的右端项b代码如下所示：
N=1000
m=10;
DIA=linspace（1,N,N）;
D=diag（DIA）;
U=orth（rand（N,N））;
A1=U'*D*U;
b1=zeros（N,1）;
fori=1:
m
b1=b1+rand
（1）*（U（i,:
））';
end
X0=zeros（N,1）;
e=10^（-6）;
tic
[S1,ek,num]=LanczosT3（A1,X0,b1,e）;
toc
对构造好的A进行计算，在计算时，取X0=zeros（N,1）;e=10-6，计算结果如下：
Arnoldi方法
GMRES方法
特征向量个数
10
20
50
100
500
750
1000
Arnoldi方法
-6.0597
-6.1419
-6.0282
-6.0041
-6.0778
-6.0089
-6.0549
迭代次数
44
87
118
132
153
155
164
运行时间
0.08051
0.17837
0.290995
0.34232
0.43846
0.45063
0.507530
GMRES
方法
-6.0244
-6.0286
-6.0567
-6.0291
-6.0311
-6.0115
-6.0510
迭代次数
43
88
112
131
148
152
160
运行时间
0.04774
0.13337
0.19433
0.25138
0.32558
0.352942
0.385624
得到以下几个结论：
1.对于基本的Arnoldi方法和GMRES方法，计算结果总是收敛的，GMRES方法收敛性较Arnoldi方法好，当m较大时，GMRES总能比Arnoldi方法提前几步收敛；
2.M较小时，计算过程中收敛误差会产生波动，m较大时，波动减小，且Arnoldi方法波动较GMRES方法剧烈；
3.收敛步数随着m值的增大而增加，但是由上表我们可以看出收敛步数的增幅逐渐减小，可以推测当m达到一定值后，收敛步数可能趋近于某一定值；
4.这个与第一次实验的Lanczos方法的结论类似。
5.
6.
7.
8.（注：
素材和资料部分来自网络，供参考。
请预览后才下载，期待你的好评与关注！
）
9.
10.
11.
升级会员 