贾哥高等数值分析第一次实验.docx
《贾哥高等数值分析第一次实验.docx》由会员分享，可在线阅读，更多相关《贾哥高等数值分析第一次实验.docx（23页珍藏版）》请在冰豆网上搜索。
贾哥高等数值分析第一次实验
高等数值分析第一次实验
第一题：
构造例子说明CG的数值形态。
当步数=阶数时CG的解如何？
当A的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时，方法的收敛性如何？
解：
用Housholder变换和对角阵构造1000阶正定对称矩阵A：
1）构造对角阵D=diag（linspace（1,1000,1000））；
2）构造Householder阵H。
取单位向量w=[1,0,0,.....0]T，I为1000阶单位矩阵，H=I–wTw。
3）构造对称正定矩阵A。
A=HTDH。
由于D是对角阵，H是对称的，所以A对称；且A与D具有相同的特征值linspace（1,1000,1000）>0，因此A对阵正定。
4）b=rand（1000,1）;取初始解x0=zeros（1000,1）;
1.计算Ax=b
利用matlab编程实现CG算法。
由于实际计算存在机器误差，因此迭代1000步后的残差不等于0，因此不能用rk=0作为停机准则，否则matlab会无休止地计算下去。
本例采用停机准则为：
迭代步数=1000步。
当D=diag（linspace（1,1000,1000））时，条件数k=1000；
当D=diag（linspace（1,100,1000））时，条件数k=100；
当D=diag（linspace（1,10,1000））时，条件数k=10；
分别计算上述三种条件数下的CG算法，得到迭代步数与残差的曲线图。
图1：
log（‖rk‖）与步数关系曲线。
横坐标是迭代步数，纵坐标是残差的对数值。
图1
如图1所示，矩阵A的条件数k越小，CG法收敛速度越快。
附matlab程序1-1：
clearall
clc
%条件数k=1000
D=diag（linspace（1,1000,1000））;
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）,'r-'）;
holdon;
k=k+1;
end
%条件数k=100
clearall
D=diag（linspace（1,100,1000））;
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）,'b-'）;
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）,'b-'）;
holdon;
k=k+1;
end
%条件数k=10
clearall
D=diag（linspace（1,10,1000））;
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）,'black-'）;
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）,'black-'）;
holdon;
k=k+1;
end
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：
clearall
clc
%条件数k=1000,特征值均匀分布
D=diag（linspace（1,1000,1000））;
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）,'r-'）;
holdon;
k=k+1;
end
%条件数k=1000，最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值
clearall
D=diag（[1,linspace（500,600,998）,1000]）;
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）,'b-'）;
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）,'b-'）;
holdon;
k=k+1;
end
title（'特征值分布对CG法收敛特性的影响'）;
xlabel（'迭代步数'）
ylabel（'残差对数log（||rk||）'）
第二题
对于同样的例子，比较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
迭代步数
时间/s
1000
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：
clearall
clc
n=1000;
D=diag（linspace（1,10,n））;%条件数分别取1000,100,10
w=eye（1,n）;
I=eye（n）;
H=I-w*w';
A=H'*D*H;%生成1000阶的对称矩阵
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（:
1）;
bita
（1）=norm（r0）;
Q（:
2）=r0/bita
（1）;%给各变量赋初始值
forj=2:
nr0=A*Q（:
j）-bita（j-1）*Q（:
j-1）;
alpha（j）=Q（:
j）'*r0;
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）;
end
fork=1:
j-1
T（k+1,k）=bita（k）;
T（k,k+1）=bita（k）;%生成三对角阵T
end
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'）;holdon;
j%迭代步数
toc
%=======CG算法======
tic
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
break;
end
beta=（R'*R）/（R0'*R0）;
p=R+beta*p;
R0=R;
end
semilogy（G,'b'）;
toc
i%迭代步数
legend（'Lanczos','CG'）
title（'Lanczos与CG算法收敛性对比,条件数k=10'）
xlabel（'迭代步数'）
ylabel（'残差对数log（||rk||）'）
第三题
当A只有m个不同特征值时，对于大的m和小的m，观察有限精度下Lanczos方法如何收敛。
解：
分别构建m=10、50、100、1000四个矩阵A，设置停机准则为：
norm（rk）图6
各个m值达到收敛时迭代步数如表2所示。
表2
m
500
800
950
990
1000
迭代步数
16
28
56
119
242
可见，相同的特征值个数越多，Lanczos收敛速度越快，而且只要m值稍小于矩阵阶数，收敛速度就明显提高，如图中m=990比m=1000的收敛速度快了一倍。
事实上，由于所构造的矩阵条件数都不坏（最大为1000），因此算法的收敛步数远小于m值。
附matlab程序3：
clearall
clc
fori=1:
5
n=1000;
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）]）;
w=eye（1,n）;
I=eye（n）;
H=I-w*w';
A=H'*D*H;%生成1000阶的对称矩阵
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（:
1）;
bita
（1）=norm（r0）;
Q（:
2）=r0/bita
（1）;%给各变量赋初始值
forj=2:
nr0=A*Q（:
j）-bita（j-1）*Q（:
j-1）;
alpha（j）=Q（:
j）'*r0;
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）;
end
fork=1:
j-1
T（k+1,k）=bita（k）;
T（k,k+1）=bita（k）;%生成三对角阵T
end
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,i）=r;
ifr/norm（b）break;
end%r的精度达到要求后停止迭代，得到最终X
end
toc
j
clearey1WXrT
end
semilogy（F（:
1）,'r'）;holdon;
semilogy（F（:
2）,'b'）;holdon;
semilogy（F（:
3）,'y'）;holdon;
semilogy（F（:
4）,'black'）;holdon;
semilogy（F（:
5）,'g'）;holdon;
legend（'m=500','m=800','m=950','m=990','m=1000'）
title（'不同特征值个数的Lanczos收敛性'）
xlabel（'迭代步数'）
ylabel（'残差对数log（||rk||）'）
第四题
取初始值近似解为零向量，右端项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
m
10
50
65
80
300
800
1000
迭代步数
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：
clearall
clc
n=1000;
D=diag（linspace（1,1000,n））;
P1=rand（n）;
[Q1,R]=qr（P1）;
A=Q1'*D*Q1;%生成1000阶的对称矩阵
[V,D1]=eig（A）;
fori=1:
7
M=[10,50,65,80,300,800,1000];
m=M（i）;
b=m*mean（V（:
1:
m）,2）;
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（:
1）;
bita
（1）=norm（r0）;
Q（:
2）=r0/bita
（1）;%给各变量赋初始值
forj=2:
nr0=A*Q（:
j）-bita（j-1）*Q（:
j-1）;
alpha（j）=Q（:
j）'*r0;
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）;
end
fork=1:
j-1
T（k+1,k）=bita（k）;
T（k,k+1）=bita（k）;%生成三对角阵T
end
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,i）=r;
ifrbreak;
end%r的精度达到要求后停止迭代，得到最终X
end
toc
j
clearey1XrT
end
semilogy（F（:
1）,'r'）;holdon;
semilogy（F（:
2）,'b'）;holdon;
semilogy（F（:
3）,'y'）;holdon;
semilogy（F（:
4）,'black'）;holdon;
semilogy（F（:
5）,'g'）;holdon;
semilogy（F（:
6）,'m'）;holdon;
semilogy（F（:
7）,'c'）;holdon;
legend（'m=10','m=50','m=65','m=80','m=300','m=800','m=1000'）
title（'b由不同特征向量个数组成时Lanczos的收敛性'）
xlabel（'迭代步数'）
ylabel（'残差对数log（||rk||）'）
第五题
构造对称不定矩阵，验证Lanczos方法的近似中断，观察收敛曲线中的峰点个数和特征值的分布关系；观测当出现峰点时，MINRES方法的收敛形态怎样。
解：
分别构建负特征值个数为m=0、10、100、500矩阵1000阶系数矩阵A，用Lanczos算法和MINRES算法求解Ax=b，设置停机准则:
norm（rk）<1e-6。
图8,9,10,11分别为m取0,10,100,500时的对比图。
图8
图9
图10
图11
通过以上实验可以看到，当负特征值个数不超过特征值总个数一半的时候，负特征值越多，峰点个数越多。
实际上，更一般的结论是，当正负特征值个数相差较多时，峰点个数较少，当正负特征值个数相当时，峰点个数最多。
Lanczos算法收敛曲线发生了非常大的波动，存在多个峰点，但相应的MINRES算法的收敛曲线却十分稳定，没有震荡。
两种算法所需的迭代次数相近（注：
m=500时均没有在1000步内达到收敛）。
附matlab程序5：
（Lanczos算法同第四题）
%======MINRES算法=====
clearall
clc
n=1000;
b=rand（n,1）;
w=eye（1,n）;
I=eye（n）;
H=I-w*w';
M=[1,10,100,500];
m=M
（1）;
D=diag（[linspace（-m,-1,m）,linspace（1,1000-m,n-m）]）;
A=H'*D*H;%生成1000阶的对称矩阵
X0=zeros（n,1）;%初始解
x=X0;
r0=b-A*X0;
R0=r0;
y=norm（r0）;%Lanczos过程
G
（1）=norm（r0）;
Q（:
1）=r0/norm（r0）;
r0=A*Q（:
1）;
alpha
（1）=Q（:
1）'*r0;
r0=r0-alpha
（1）*Q（:
1）;
bita
（1）=norm（r0）;
Q（:
2）=r0/bita
（1）;%给各变量赋初始值
forj=2:
500r0=A*Q（:
j）-bita（j-1）*Q（:
j-1）;
alpha（j）=Q（:
j）'*r0;
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）;
end
fork=1:
j-1
T（k+1,k）=bita（k）;
T（k,k+1）=bita（k）;%生成三对角阵T
end
e
（1）=1;
e（2:
j）=0;
%对T进行QR分解
[Qk,Rk]=qr（T）;
gk=Qk'*norm（r0）*e';
yk=Rk\gk;
xk=Q（:
1:
j）*yk;
rk=abs（gk（end））;
G（j）=rk;
ifrk<1e-6;
升级会员 