清华大学高等数值分析第一次实验.docx
- 文档编号:3927605
- 上传时间:2022-11-26
- 格式:DOCX
- 页数:28
- 大小:541.15KB
清华大学高等数值分析第一次实验.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));
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 清华大学 高等 数值 分析 第一次 实验