贾哥高等数值分析第一次实验.docx
- 文档编号:8639641
- 上传时间:2023-02-01
- 格式:DOCX
- 页数:23
- 大小:263.24KB
贾哥高等数值分析第一次实验.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; ifr break; 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;
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 高等 数值 分析 第一次 实验
![提示](https://static.bdocx.com/images/bang_tan.gif)