贾哥高等数值分析报告第一次实验Word格式.docx
《贾哥高等数值分析报告第一次实验Word格式.docx》由会员分享，可在线阅读，更多相关《贾哥高等数值分析报告第一次实验Word格式.docx（23页珍藏版）》请在冰豆网上搜索。
w=eye（1,1000）;
I=eye（1000）;
H=I-w*w'
;
A=H'
*D*H;
%生成1000阶的对称矩阵
b=rand（1000,1）;
x=zeros（1000,1）;
%初始解
r=b-A*x;
%初始残量
p=r;
%初始搜索方向
k=0;
semilogy（0,norm（r）,'
r-'
）;
holdon;
whilek<
1000
alpha=r'
*p/（p'
*A*p）;
x=x+alpha*p;
rold=r;
r=rold-alpha*A*p;
beta=r'
*r/（rold'
*rold）;
p=r+beta*p;
semilogy（k,norm（r）,'
holdon;
k=k+1;
end
%条件数k=100
D=diag（linspace（1,100,1000））;
b-'
%条件数k=10
D=diag（linspace（1,10,1000））;
black-'
title（'
条件数的大小对CG法收敛特性的影响'
xlabel（'
迭代步数'
）
ylabel（'
残差对数log（||rk||）'
2.构造特殊特征值分布
构造对称正定矩阵A1和A2。
D1=diag（linspace（1,1000,1000））时，条件数k=1000，特征值均匀分布；
D2=diag（[1,linspace（500,600,998）,1000]）时，条件数仍为k=1000，最大特征值1000远大于第二个最大特征值600，最小特征值1远小于第二个最小特征值500。
图2：
矩阵特征值分布对CG算法收敛性的影响
图2
如图2所示，A1和A2的条件数均为1000，但A2的收敛速度远高于A1。
这是因为，在CG算法中，系数矩阵的中间特征值分布对CG的收敛速度有巨大的影响。
经过几步后，CG的收敛因子将是：
=0.046
而非：
=0.939
因此，A2矩阵的收敛速度快得多。
附matlab程序1-2：
%条件数k=1000,特征值均匀分布
%条件数k=1000，最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值
D=diag（[1,linspace（500,600,998）,1000]）;
特征值分布对CG法收敛特性的影响'
第二题
对于同样的例子，比较CG和Lanczos的计算结果
采用和第一题相同的构造方法，构造三个1000阶正定对称矩阵，使条件数k分别为：
1000,100,10。
分别采用CG和Lanczos方法计算Ax=b，且都设置停机准则为：
norm（rk）<
1e-8。
图3，图4，图5分别为3种条件数下的CG法与Lanczos法比较。
图3
图4
图5
两种算法的收敛情况见表1。
表1
对比
Lanczos
CG
条件数
迭代步数
时间/s
243
0.82
267
0.41
100
124
0.43
134
0.21
10
41
0.24
43
0.07
可以看到，Lanczos的迭代步数比CG小，说明Lanczos收敛更快，但是其运算时间却明显大于CG，因此总体上CG的收敛效率要比Lanczos高。
因此在计算对称正定问题时，优先选用CG算法。
附matlab程序2：
n=1000;
D=diag（linspace（1,10,n））;
%条件数分别取1000,100,10
w=eye（1,n）;
I=eye（n）;
b=rand（n,1）;
X0=zeros（n,1）;
x=X0;
r0=b-A*X0;
R0=r0;
%======Lanczos解法========
tic
y=norm（r0）;
F
（1）=norm（r0）;
Q（:
1）=r0/norm（r0）;
r0=A*Q（:
1）;
alpha
（1）=Q（:
1）'
*r0;
r0=r0-alpha
（1）*Q（:
bita
（1）=norm（r0）;
2）=r0/bita
（1）;
%给各变量赋初始值
forj=2:
nr0=A*Q（:
j）-bita（j-1）*Q（:
j-1）;
alpha（j）=Q（:
j）'
r0=r0-alpha（j）*Q（:
j）;
ifnorm（r0）>
1
bita（j）=norm（r0）;
Q（:
j+1）=r0/bita（j）;
end
fork=1:
j
T（k,k）=alpha（k）;
j-1
T（k+1,k）=bita（k）;
T（k,k+1）=bita（k）;
%生成三对角阵T
e
（1）=1;
e（2:
j）=0;
y1=T\（y*e）'
W=Q（:
1:
j）;
X=X0+W*y1;
r=norm（b-A*X）;
%求解第k步生成的X及r
F（j）=r;
ifr/norm（b）<
1e-8
break;
end%r的精度达到要求后停止迭代，得到最终X
end
semilogy（F,'
r'
j%迭代步数
toc
%=======CG算法======
p=R0;
fori=1:
n
R=R0;
a=（R'
*R）/（p'
*（A*p））;
x=x+a*p;
R=R-a*（A*p）;
G（i）=norm（R）;
ifG（i）/norm（b）<
=1e-8
beta=（R'
*R）/（R0'
*R0）;
p=R+beta*p;
R0=R;
semilogy（G,'
b'
i%迭代步数
legend（'
Lanczos'
'
CG'
Lanczos与CG算法收敛性对比,条件数k=10'
第三题
当A只有m个不同特征值时，对于大的m和小的m，观察有限精度下Lanczos方法如何收敛。
分别构建m=10、50、100、1000四个矩阵A，设置停机准则为：
sqrt（eps），分别计算Ax=b的解，如图6所示。
图6
各个m值达到收敛时迭代步数如表2所示。
表2
m
500
800
950
990
16
28
56
119
242
可见，相同的特征值个数越多，Lanczos收敛速度越快，而且只要m值稍小于矩阵阶数，收敛速度就明显提高，如图中m=990比m=1000的收敛速度快了一倍。
事实上，由于所构造的矩阵条件数都不坏（最大为1000），因此算法的收敛步数远小于m值。
附matlab程序3：
5
M=[500,800,950,990,1000];
m=M（i）;
%m=500、800、950、990、1000
D=diag（[（1001-m）*ones（1,n-m）,linspace（1001-m,1000,m）]）;
F（j,i）=r;
sqrt（eps）
clearey1WXrT
semilogy（F（:
1）,'
2）,'
3）,'
y'
4）,'
black'
5）,'
g'
m=500'
m=800'
m=950'
m=990'
m=1000'
不同特征值个数的Lanczos收敛性'
第四题
取初始值近似解为零向量，右端项b仅由A的m个不同特征向量的线性组合表示时，Lanczos方法的收敛性如何？
数值计算中方法的收敛性和m的大小关系如何？
对任意一个随机的矩阵P进行QR分解，可以得到一个正交阵Q，再令A=QT*D*Q得到系数矩阵A。
则b可以由Q的前m个列向量相加得到。
停机标准：
其matlab程序与第三题基本相同，只是前面构造A和b的方式略有不同，附matlab程序5中只写出了不同的部分。
b由m个不同特征向量组成时的Lanczos收敛情况如图7所示。
图7
m取不同值时的收敛步数见表3。
表3
50
65
80
300
6
8
68
129
169
188
192
可见，随着m的增大，迭代步数逐渐增加。
理论上，当b仅由A的m个不同特征向量的线性组合表示时，Lanczos方法必至多m步收敛。
但实际上，当m较大时，由于精度等方面的限制，m个特征向量并不都线性无关，所以，当m较大时，迭代步数可能比m小的多，比如本例中m>
300后，迭代步数远小于m。
而同样由于计算精度限制，在m较小时，迭代步数可能溢出m步，比如本例m=65和80的情形。
另外，m值大于一定值后，随着m的增加，收敛步数逐渐趋于定值，可能的原因是：
随着m的增大，特征向量线性相关性逐渐增强。
附matlab程序4：
D=diag（linspace（1,1000,n））;
P1=rand（n）;
[Q1,R]=qr（P1）;
A=Q1'
*D*Q1;
[V,D1]=eig（A）;
7
M=[10,50,65,80,300,800,1000];
b=m*mean（V（:
m）,2）;
ifr<
clearey1XrT
6）,'
m'
7）,'
c'
m=10'
m=50'
m=65'
m=80'
m=300'
b由不同特征向量个数组成时Lanczos的收敛性'
第五题
构造对称不定矩阵，验证Lanczos方法的近似中断，观察收敛曲线中的峰点个数和特征值的分布关系；
观测当出现峰点时，MINRES方法的收敛形态怎样。
分别构建负特征值个数为m=0、10、100、500矩阵1000阶系数矩阵A，用Lanczos算法和MINRES算法求解Ax=b，设置停机准则:
1e-6。
图8,9,10,11分别为m取0,10,100,500时的对比图。
图8
图9
图10
图11
通过以上实验可以看到，当负特征值个数不超过特征值总个数一半的时候，负特征值越多，峰点个数越多。
实际上，更一般的结论是，当正负特征值个数相差较多时，峰点个数较少，当正负特征值个数相当时，峰点个数最多。
Lanczos算法收敛曲线发生了非常大的波动，存在多个峰点，但相应的MINRES算法的收敛曲线却十分稳定，没有震荡。
两种算法所需的迭代次数相近（注：
m=500时均没有在1000步内达到收敛）。
附matlab程序5：
（Lanczos算法同第四题）
%======MINRES算法=====
M=[1,10,100,500];
m=M
（1）;
D=diag（[linspace（-m,-1,m）,linspace（1,1000-m,n-m）]）;
A=H'
X0=zeros（n,1）;
x=X0;
r0=b-A*X0;
R0=r0;
y=norm（r0）;
%Lanczos过程
G
（1）=norm（r0）;
r0=A*Q（:
alpha
（1）=Q（:
r0=r0-alpha
（1）*Q（:
bita
（1）=norm（r0）;
500r0=A*Q（:
%对T进行QR分解
[Qk,Rk]=qr（T）;
gk=Qk'
*norm（r0）*e'
yk=Rk\gk;
xk=Q（:
j）*yk;
rk=abs（gk（end））;
G（j）=rk;
ifrk<
1e-
升级会员 