清华大学高等数值分析第一次实验.docx
《清华大学高等数值分析第一次实验.docx》由会员分享，可在线阅读，更多相关《清华大学高等数值分析第一次实验.docx（28页珍藏版）》请在冰豆网上搜索。
清华大学高等数值分析第一次实验
材料科学与工程学院
高等数值分析实验一
姓名：
薛康乐
学号：
2013211320
2013/12/6
高等数值分析第一次实验
班级：
材硕13班学号：
**********姓名：
薛康乐
一、构造例子说明CG的数值性态。
当步数=阶数时CG的解如何？
当A的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时方法的收敛性如何？
1、构造步数=阶数的CG算法例子
（1）首先随机生成n阶矩阵A1，若A1不满秩，则可由此构造对称正定矩阵A=A1’*A1，设定精确解Xe为元素全部为1的n维列向量，则b=A*Xe。
（2）进行CG算法求解方程组，分别取n=1000,2000,3000,4000进行计算，获得不同阶数下的运行时间、相对残差和相对误差。
由于方程组是随机给出的，所以同时计算不同方程组条件下的条件数。
（3）计算结果如下所示：
方程阶数n
算法时间（s）
相对残差
相对误差
方程组条件数
1000
6.40E+00
3.32E-11
5.39E-04
1.03E+09
2000
54.440142
2.43E-11
1.23E-04
1.22E+10
3000
168.762032
5.29E-13
4.31E-04
3.03E+10
4000
413.927374
1.07E-12
9.70E-05
4.86E+11
1000阶方程残差随迭代次数变化关系（左:
rk，右:
lg（rk））
2000阶方程残差随迭代次数变化关系（左:
rk，右:
lg（rk））
3000阶方程残差随迭代次数变化关系（左:
rk，右:
lg（rk））
4000阶方程残差随迭代次数变化关系（左:
rk，右:
lg（rk））
CG算法在步数小于阶数是即可取得稳定解。
对数坐标下的收敛曲线不光滑，震荡严重，残差的2范数没有最优性。
但是方法的相对残差大体上成下降趋势，最终仍能够达到收敛。
但是随着阶数的增加，总体残差是变小的，但是变小的趋势又是减缓的。
Matlab程序如下：
>>clearall
n=1000;
A1=zeros（n,n）;
Xe=ones（n,1）;
whilerank（A1）~=n;
A1=randi（[0,10],n,n）;
end
A=A1'*A1;
b=A*Xe;
tic
X0=zeros（n,1）;
ro=b-A*X0;
p1=ro;
r1=ro;
fork=1:
n
ro=r1;
alpha=（ro'*ro）/（p1'*（A*p1））;
X0=X0+alpha*p1;
r1=ro-alpha*（A*p1）;
e（k）=sqrt（r1'*r1）/（sqrt（b'*b）+norm（A1,1）*sqrt（X0'*X0））;
bita=（r1'*r1）/（ro'*ro）;
p1=r1+bita*p1;
end;
toc
e（n）
sqrt（（Xe-X0）'*（Xe-X0））/sqrt（Xe'*Xe）
cond（A）
subplot（1,2,1）;plot（[1:
k],e）;
subplot（1,2,2）;plot（[1:
k],log（e））;
2、当矩阵为题设中的特殊矩阵时
（1）首先随机生成n阶矩阵B，并保证B不满秩。
构造对角矩阵A1（最大特征值远大于第二大特征值）和A2（最小特征值远小于第二小特征值），则A可以由b1=B’*B,b2=b1’*A1*b1,A=b2’*A2*b2。
设定精确解Xe为元素全部为1的n维列向量，则b=A*Xe。
（2）编程计算
利用CG算法进行求解，分别另n=1000,2000,3000,4000。
获得不同阶数下的运行时间、相对残差和相对误差。
由于方程组是随机给出的，所以同时计算不同方程组条件下的条件数。
（3）计算结果如下所示：
方程阶数n
算法时间（s）
相对残差
相对误差
方程组条件数
1000
6.497919
6.93E-12
0.0043
6.18E+18
2000
49.407932
1.10E-14
0.0032
2.68E+19
3000
178.003284
2.11E-14
0.0073
2.82E+19
4000
369.241012
4.43E-11
0.0257
2.16E+19
1000阶方程残差随迭代次数变化关系（左:
rk，右:
lg（rk））
2000阶方程残差随迭代次数变化关系（左:
rk，右:
lg（rk））
3000阶方程残差随迭代次数变化关系（左:
rk，右:
lg（rk））
结论：
随着方程阶数的增高，收敛速度变慢，计算时间与阶数的三次方成正比。
收敛过程中会出现残差严重震荡的情况。
而且随着方程阶数的增加，会出现残差随着迭代进行整体增大的情况（如上图所示）。
并且与特征值没有特殊分布的情况相比，特征值存在特殊分布的方程组的残差离散程度更大。
Matlab程序如下所示：
>>clearall
n=4000;
a=1:
2:
8000;
a（n）=1e7;
a
（1）=1e-5;
B1=zeros（n,n）;
Xe=ones（n,1）;
whilerank（B1）~=n
B1=randi（[0,10],n,n）;
end;
A1=diag（a）;
b1=B1'*B1;
A=b1'*A1*b1;
b=A*Xe;
tic
X0=zeros（n,1）;
ro=b-A*X0;
p1=ro;
r1=ro;
fork=1:
n
ro=r1;
alpha=（ro'*ro）/（p1'*（A*p1））;
X0=X0+alpha*p1;
r1=ro-alpha*（A*p1）;
e（k）=sqrt（r1'*r1）/（sqrt（b'*b）+norm（A1,1）*sqrt（X0'*X0））;
bita=（r1'*r1）/（ro'*ro）;
p1=r1+bita*p1;
end;
toc
e（n）
sqrt（（Xe-X0）'*（Xe-X0））/sqrt（Xe'*Xe）
subplot（1,2,1）;plot（[1:
k],e）;
subplot（1,2,2）;plot（[1:
k],log（e））;
cond（A）
二、对于同样的例子，比较CG和Lanczos的计算结果
（1）构造方程组，首先创建对称正定矩阵A，随机生成一个满秩的n维矩阵M，则A为对称正定矩阵。
设定精确解Xe为元素全部为1的n维列向量，则b=A*Xe。
（2）编程计算。
令n=2000，分别用CG和Lanczos算法进行计算，比较其收敛步数、相对残差和相对误差。
（3）结果及分析
试验方法
Lanczos法
CG法
收敛步数
144
146
绝对误差
-18.3095
0.0135
相对误差
9.7698e-09
9.6989e-09
lanczos法的残差（左:
rk，右:
lg（rk））
lanczos法残差下降很快，但是通常几部后qj的数值正交性就会失去，最终的相对残差也会比CG的大，残差的震荡程度比CG法要小。
收敛步数与CG法相似。
CG法残差（左:
rk，右:
lg（rk））
而相同阶数的CG法，就收敛速度而言，开始确实比Lanczos法要慢，但是整体收敛速度比Lanczos好。
起始阶段，CG法的残差震荡很严重，但最后获得的相对残差，比Lanczos法获得的相对残差要小。
Matlab程序如下所示：
clearall
n=2000;
A1=zeros（n,n）;
Xe=ones（n,1）;
whilerank（A1）~=n;
A1=randi（[0,10],n,n）;
end
A=A1'*A1;
b=A*Xe;
X0=rand（n,1）;
r0=b-A*X0;
x=X0;
X=X0;
R0=r0;
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:
n
r0=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）;
end
e
（1）=1;
e（2:
j）=0;
y1=T\（y*e）';
r=bita（j）*abs（y1（end））;
F（j）=r;
ifr/norm（b）<1e-8
break;
end
end
figure;
subplot（1,2,1）;plot（F）;
subplot（1,2,2）;plot（log（F））;
j
norm（X）-norm（Xe）
r/norm（b）
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
be=（R'*R）/（R0'*R0）;
p=R+be*p;
R0=R;
end
re=sqrt（（Xe-x）'*（Xe-x））/sqrt（Xe'*Xe）;
figure;
subplot（1,2,1）;plot（G）;
subplot（1,2,2）;plot（log（G））;
norm（x）-norm（Xe）
G（i）/norm（b）
三、当A只有m个不同特征值时，对于大的m和小的m，观察有限精度下Lanczos方法如何收敛。
（1）先构造不同特征值个数分别为1,10,100,500,1000,2000的n阶矩阵M，通过对某n维随机矩阵P进行QR分解得到正交阵Q，则由谱分解原理可知，n阶矩阵A=Q’*M*Q就是满足条件的对称正定矩阵。
设精确解Xe为元素全部为1的n维列向量，则b=A*Xe。
（2）编程计算。
（3）结果及分析
不同的m
1
10
100
500
1000
2000
迭代步数
2
10
52
110
155
204
相对残差
0
-41.2997
-33.8495
-21.7892
-12.6551
-1.9448e-10
相对误差
3.5197e-30
2.5151e-15
8.4156e-09
9.1810e-09
8.8403e-09
9.6940e-09
m=1时的收敛曲线（左:
rk，右:
lg（rk））
m=10时的收敛曲线（左:
rk，右:
lg（rk））
m=100时的收敛曲线（左:
rk，右:
lg（rk））
m=500时的收敛曲线（左:
rk，右:
lg（rk））
m=1000时的收敛曲线（左:
rk，右:
lg（rk））
m=2000时的收敛曲线（左:
rk，右:
lg（rk））
从上述结果中可以看出，m值越大，收敛越慢，当矩阵的特征值越密集时收敛越快。
从相对收敛速度而言，当m大时，虽然收敛时的迭代步数较大，但是这个步数与m相比还是较小的；当m小时，虽然收敛时的迭代步数较小，但是这个步数与m相比还是较大的。
matlab程序如下所示：
n=2000;
m=1;
v=ones（1,2000）;
（当m不等于1时，v由下式给出：
v=horzcat（[（n/m）:
（n/m）:
n],zeros（1,（n-m））;）
D=diag（v）;
P=rand（n）;
[Q,R]=qr（P）;
A=Q'*D*Q;
Xe=ones（n,1）;
b=A*Xe;
X0=zeros（n,1）;
r0=b-A*X0;
x=X0;
R0=r0;
y=norm（r0）;
F
（1）=norm（r0）;
S（:
1）=r0/norm（r0）;
r0=A*S（:
1）;
alpha
（1）=S（:
1）'*r0;
r0=r0-alpha
（1）*S（:
1）;
bita
（1）=norm（r0）;
S（:
2）=r0/bita
（1）;
forj=2:
n
r0=A*S（:
j）-bita（j-1）*S（:
j-1）;
alpha（j）=S（:
j）'*r0;
r0=r0-alpha（j）*S（:
j）;
ifnorm（r0）>0
bita（j）=norm（r0）;
S（:
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）;
end
e
（1）=1;
e（2:
j）=0;
y1=T\（y*e）';
W=S（:
1:
j）;
X=X0+W*y1;
r=bita（j）*abs（y1（end））;
F（j）=r;
ifr/norm（b）<1e-8
break;
end
end
figure;
subplot（1,2,1）;plot（F/norm（b））;
subplot（1,2,2）;plot（log（F/norm（b）））;
j
norm（X）-norm（Xe）
r/norm（b）
四、取初始近似解为零向量，右端项b仅由A的m个不同个特征向量的线性组合表示时，Lanczos方法的收敛性如何？
数值计算中方法的收敛性和m的大小关系如何：
（1）构造算例：
先构造n阶对角矩阵M，通过对某n维随机矩阵P进行QR分解之后得到正交阵Q，则由谱分解原理可知n阶矩阵A=Q’*M*Q就是满足条件的对称正定矩阵；通过对A进行特征值分解得到其特征向量，取其中m个进行线性组合，则右端项b=A*Xe。
设定精确解为元素全部为1的n维列向量。
（2）编程计算：
取n=2000，m=1、10、100、500、1000、1500、2000分别利用Lanczos方法进行计算，并对各情况下的收敛步数、相对残差和相对误差进行分析比较。
（3）计算结果及分析
m
1
10
100
500
1000
1500
2000
迭代步数
2
3
139
192
203
209
213
相对残差
-2.8311e-14
3.1086e-15
-1.4837e-10
-1.3723e-10
-1.5347e-10
-1.5935e-10
-1.7764e-10
相对误差
1.4514e-09
2.2118e-09
9.8056e-09
9.0426e-09
9.0195e-09
8.5207e-09
9.1195e-09
m=1时的收敛曲线（左:
rk，右:
lg（rk））
m=10时的收敛曲线（左:
rk，右:
lg（rk））
m=100时的收敛曲线（左:
rk，右:
lg（rk））
m=1000时的收敛曲线（左:
rk，右:
lg（rk））
m=1500时的收敛曲线（左:
rk，右:
lg（rk））
m=2000时的收敛曲线（左:
rk，右:
lg（rk））
收敛次数与m的关系
从上述图表可以看出：
收敛特性与特征向量相关。
当m增大时，迭代次数随之增大，但是收敛步数增大的程度小于m增大的程度，当m达到一定值后，收敛步数可能趋于某一定值。
随着m的增加，相对残差也呈现增加的趋势，但在确定精度要求的前提下，收敛效果随着m的增加而增加。
matlab程序如下所示：
n=2000;
m=2000;%根据需要改变m的值
v=[1:
1:
n];
D=diag（v）;
P=rand（n）;
[Q,R]=qr（P）;
A=Q'*D*Q;
[V,D]=eig（A）;
Xe=zeros（n,1）;
fori=1:
m
Xe=Xe+V（:
i）;
end
b=A*Xe;
X0=zeros（n,1）;
r0=b-A*X0;
x=X0;
R0=r0;
y=norm（r0）;
F
（1）=norm（r0）;
S（:
1）=r0/norm（r0）;
r0=A*S（:
1）;
alpha
（1）=S（:
1）'*r0;
r0=r0-alpha
（1）*S（:
1）;
bita
（1）=norm（r0）;
S（:
2）=r0/bita
（1）;
forj=2:
n
r0=A*S（:
j）-bita（j-1）*S（:
j-1）;
alpha（j）=S（:
j）'*r0;
r0=r0-alpha（j）*S（:
j）;
ifnorm（r0）>0
bita（j）=norm（r0）;
S（:
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）;
end
e
（1）=1;
e（2:
j）=0;
y1=T\（y*e）';
W=S（:
1:
j）;
X=X0+W*y1;
r=bita（j）*abs（y1（end））;
F（j）=r;
ifr/norm（b）<1e-8
break;
end
end
figure;
subplot（1,2,1）;plot（F/norm（b））;
subplot（1,2,2）;plot（log（F/norm（b）））;
j
norm（X）-norm（Xe）
r/norm（b）
五、构造对称不定的矩阵，验证Lanczos方法的近似中断，观察收敛曲线中的峰点个数和特征值的分布关系，观察当出现峰点时，MINRES方法的收敛性态怎样。
（1）构造算例：
首先构造6个n阶对角阵D，其对角线上分别有m=1,10,50,100,200,500个负值（既有相应个数的负特征根），其余特征根分两种情况。
A、有m个为正数，其余的都为零；B、全部为正数。
通过对某n维随机矩阵QR分解获得正交阵Q，则由谱分解原理可知，n阶矩阵A=Q’*D*Q即满足条件的对称矩阵。
设定精确解Xe为元素全部为1的n维列向量。
则b=A*Xe。
（2）编程计算：
取n=2000，分别利用Lanczos方法进行计算。
（3）计算结果及分析
A、有m个为正数，其余的都为零
负特征值个数m
1
20
50
100
200
500
Lanczos收敛步数
2
44
126
323
638
1559
相对残差
4.3215e-15
3.8191e-09
2.4956e-09
3.0200e-09
8.0355e-09
9.7086e-09
Minres收敛步数
2
44
125
320
613
1540
相对残差
2.6e-15
5.5e-09
9.8e-09
9.9e-09
9.5e-09
9.8e-09
m=1的收敛曲线（左：
Lanczos，右：
MINRES）
m=20的收敛曲线（左：
Lanczos，右：
MINRES）
m=50的收敛曲线（左：
Lanczos，右：
MINRES）
m=100的收敛曲线（左：
Lanczos，右：
MINRES）
m=200的收敛曲线（左：
Lanczos，右：
MINRES）
m=500的收敛曲线（左：
Lanczos，右：
MINRES）
当特征值中含有0时，迭代过程很快中断，并且负特征值越多，Lanczos方法的峰点越多，震荡越严重。
MINRES方法与Lanczos方法收敛的步数相近，但是收敛曲线非常稳定，没有震荡。
在实验中，所选用的m值都在2000步之内获得了收敛。
MINRES方法的收敛性态比Lanczos方法要好。
为了进行对比，又进行了一组实验，其情况B。
具体如下所示：
B、有m个为负数，其余的都为正数
负特征值个数m
1
20
50
100
200
500
Lanczos收敛步数
228
1535
2000（未收敛）
2000（未收敛）
2000（未收敛）
2000（未收敛）
相对残差
8.7500e-09
7.8854e-09
2.3168e-06
9.2077e-06
2.8337e-04
5.7899e-04
Minres收敛步数
214
1304
2000（未收敛）
2000（未收敛）
2000（未收敛）
2000（未收敛）
相对残差
9.7e-09
9.8e-09
5.1e-08
3.9e-07
2.5e-06
1.9e-05
m=1的收敛曲线（左：
Lanczos，右：
MINRES）
m=20的收敛曲线（左：
Lanczos，右：
MINRES）
m=50的收敛曲线（左：
Lanczos，右：
MINRES）
m=100的收敛曲线（左：
Lanczos，右：
MINRES）
m=200的收敛曲线（左：
Lanczos，右：
MINRES）
m=500的收敛曲线（左：
Lanczos，右：
MINRES）
当A的特征值里面不存在0值时，收敛速度明显放慢。
而且当m大到一定程度之后，两种方法都无法在要求的精度范围内收敛。
Lanczos方法依旧存在很严重的震荡，而且震荡程度随着m的增加而增加。
MINRES方法收敛曲线平滑，具有较好的收敛性态。
Matlab代码如下：
clearall;
n=2000;
v=rand（1,n）;
M=500;
fori=1:
M
v（i）=-100*rand
（1）;
v（M+i+2）=200*rand
（1）;
end
D=diag（v）;
P=rand（n）;
[Q,R]=qr（P）;
A=Q'*D*Q;
[V,D]=eig（A）;
Xe=ones（n,1）;
b=A*Xe;
X0=zeros（n,1）;
r0=b-A*X0;
R0=r0;
y=norm（r0）;
F
（1）=norm（r0）;
S（:
1）=r0/norm（r0）;
r0=A*S（:
1）;
alpha
（1）=S（:
1）'*r0;
r0=r0-alpha
（1）*S（:
1）;
bita
（1）=norm（r0）;
S（:
2）=r0/bita
（1）;
forj=2:
n
r0=A*S（:
j）-bita（j-1）*S（:
j-1）;
alpha（j）=S（:
j）'*r0;
r0=r0-alpha（j）*S（:
j）;
ifnorm（r0）>0
bita（j）=norm（r0）;
S（:
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）;
end
e
（1）=1;
e（2:
j）=0;
y1=T\（y*e）';
W=S（:
1:
j）;
X=X0+W*y1;
r=bita（j）*abs（y1（end））;
升级会员 