贾哥高等数值分析第一次实验0001

1、高等数值分析第一次实验第一题：构造例子说明 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

2、 对称；且 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( lin

3、space(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：clear allclc%条件数 k=1000D=diag(linspace(1,1000,1000);w=eye(1,1000);I=eye(1000);H=I-w*w&

4、#39;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-' );hold on;while k<1000alpha = 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), '

5、r-' );hold on ;k=k+1;end%条件数 k=100clear allD=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-' );hold on;while k<1000alpha = r'*p/(p'*A*p);x =

6、x+alpha*p;rold = r;'b-' );r = rold-alpha*A*p; beta = r'*r/(rold'*rold); p = r+beta*p; semilogy(k,norm(r), hold on ;k=k+1;end %条件数 k=10 clear allD=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; %

7、初始残量 p=r; %初始搜索方向k=0; semilogy(0,norm(r), hold on;while k<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), hold on ; k=k+1;end title( xlabel( ylabel('black-' );'black-' );条件数的大小对

8、9; 迭代步数 ' )' 残差对数 log(|rk|)'CG 法收敛特性的影响 ' )2.构造特殊特征值分布构造对称正定矩阵 A1 和 A2 。D1=diag( linspace(1, 1000, 1000) ) 时，条件数 k=1000 ，特征值均匀分布；1000 远大D2=diag(1,linspace(500,600,998),1000) 时，条件数仍为 k=1000 ，最大特征值 于第二个最大特征值 600，最小特征值 1 远小于第二个最小特征值 500。 图 2：矩阵特征值分布对 CG 算法收敛性的影响图2如图 2 所示， A1 和 A2 的条件数均为

9、 1000，但 A2 的收敛速度远高于 A1 。这是因为，在 CG 算法中，系数矩阵的中间特征值分布对 CG 的收敛速度有巨大的影响。经过几步后， CG 的收敛因子将是：而非：1=0.9391因此， A2 矩阵的收敛速度快得多。附 matlab 程序 1-2：clear all 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

10、-A*x; %初始残量p=r; %初始搜索方向k=0;semilogy(0,norm(r), 'r-' );hold on;while k<1000alpha = 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-' );hold on ;k=k+1;end%条件数 k=1000 ，最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特

11、征值 clear allD=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-' );hold on;while k<1000alpha = r'*p/(p'*A*p);x = x+alpha*p;rold = r;r = rold

12、-alpha*A*p;beta = r'*r/(rold'*rold);p = r+beta*p;semilogy(k,norm(r), 'b-' );holdon ;k=k+1;endtitle( '特征值分布对 CG 法收敛特性的影响xlabel(' 迭代步数 ' )ylabel(' 残差对数 log(|rk|)')第二题 对于同样的例子，比较 CG 和 Lanczos 的计算结果 解：采用和第一题相同的构造方法，构造三个 1000 阶正定对称矩阵，使条件数 k 分别为： 1000,100,10。分别采用 CG 和 L

13、anczos 方法计算 Ax=b ，且都设置停机准则为： norm(rk)<1e-8 。 图 3，图 4，图 5 分别为 3 种条件数下的 CG 法与 Lanczos 法比较。图3图4图5两种算法的收敛情况见表 1。 表1对比LanczosCG条件数迭代步数时间 /s迭代步数时间 /s10002430.822670.411001240.431340.2110410.24430.07可以看到， Lanczos 的迭代步数比 CG 小，说明 Lanczos 收敛更快，但是其运算时间却 明显大于 CG ，因此总体上 CG 的收敛效率要比 Lanczos 高。因此在计算对称正定问题时， 优先选用

14、 CG 算法。附 matlab 程序 2： clear all clc n=1000;D=diag(linspace(1,10,n);%条件数分别取 1000,100,10w=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 解法 = ticy=norm(r0);F(1)=norm(r0);Q(:,1)=r0/norm(r0);r0=A*Q(:,1);alpha(1)=Q(:,1)'*r0;r0=r

15、0-alpha(1)*Q(:,1);bita(1)=norm(r0);Q(:,2)=r0/bita(1);%给各变量赋初始值for j=2:n r0=A*Q(:,j)-bita(j-1)*Q(:,j-1);alpha(j)=Q(:,j)'*r0;r0=r0-alpha(j)*Q(:,j);if norm(r0)>1bita(j)=norm(r0);Q(:,j+1)=r0/bita(j);endfor k=1:jT(k,k)=alpha(k);endfor k=1:j-1T(k+1,k)=bita(k);T(k,k+1)=bita(k);%生成三对角阵 Tende(1)=1;e(2

16、:j)=0;y1=T(y*e)'W=Q(:,1:j);X=X0+W*y1;r=norm(b-A*X);%求解第 k 步生成的 X 及 rF(j)=r;if r/norm(b)<1e-8break ;end %r 的精度达到要求后停止迭代，得到最终 X endsemilogy(F, 'r' );hold on;j %迭代步数toc%=CG 算法 = ticp=R0;for i=1:nR=R0;a=(R'*R)/(p'*(A*p);x=x+a*p;R=R-a*(A*p);G(i)=norm(R);if G(i)/norm(b)<=1e-8 bre

17、ak ;end beta=(R'*R)/(R0'*R0); p=R+beta*p; R0=R;end semilogy(G, 'b' ); toci%迭代步数legend( 'Lanczos' , 'CG' ) title( 'Lanczos 与 CG 算法收敛性对比 , 条件数 k=10' ) xlabel( ' 迭代步数 ' )ylabel( ' 残差对数 log(|rk|)')第三题Lanczos 方法如何收norm(rk)<sqrt(eps) ，当A只有 m个不同特征值

18、时，对于大的 m和小的 m，观察有限精度下 敛。解：分别构建 m = 10 、 50、100、1000 四个矩阵 A，设置停机准则为： 分别计算 Ax = b 的解，如图 6 所示。图6各个 m 值达到收敛时迭代步数如表 2 所示。 表2m5008009509901000迭代步数162856119242可见，相同的特征值个数越多， Lanczos 收敛速度越快， 而且只要 m 值稍小于矩阵阶数， 收敛速度就明显提高，如图中 m=990 比 m=1000 的收敛速度快了一倍。事实上，由于所构造的矩阵条件数都不坏(最大为1000)，因此算法的收敛步数远小于m 值。附 matlab 程序 3：cle

19、ar allclcfor i=1:5n=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 解法 =ticy=norm(r0);F(1)=norm(r0);Q(:,

20、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);%给各变量赋初始值for j=2:n r0=A*Q(:,j)-bita(j-1)*Q(:,j-1);alpha(j)=Q(:,j)'*r0;r0=r0-alpha(j)*Q(:,j);if norm(r0)>1bita(j)=norm(r0);Q(:,j+1)=r0/bita(j);endfor k=1:jT(k,k)=alpha(k);endfor k=1:j-1

21、T(k+1,k)=bita(k);T(k,k+1)=bita(k);end%生成三对角阵 Te(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 及 rF(j,i)=r;if r/norm(b)<sqrt(eps)breakend %r 的精度达到要求后停止迭代，得到最终endtocclear e y1 W X r T endsemilogy(F(:,1),'r' );holdonsemilogy(F(:,2),'b' );holdonsemilog

22、y(F(:,3),'y' );holdonsemilogy(F(:,4),'black');holdsemilogy(F(:,5),'g' );holdonon;legend( 'm=500', 'm=800' , 'm=950', 'm=990', 'm=1000' )title（ ' 不同特征值个数的 Lanczos 收敛性 ' ） xlabel（ ' 迭代步数 ' ）ylabel( ' 残差对数 log(|rk|)'

23、; )第四题取初始值近似解为零向量， 右端项 b 仅由 A 的 m 个不同特征向量的线性组合表示时， Lanczos 方法的收敛性如何？数值计算中方法的收敛性和 m 的大小关系如何？解：对任意一个随机的矩阵 P进行 QR分解，可以得到一个正交阵 Q，再令 A=Q T*D*Q 得到系 数矩阵 A。则 b 可以由 Q 的前 m 个列向量相加得到。停机标准：其 matlab 程序与第三题基本相同，只是前面构造A 和 b 的方式略有不同，附 matlab 程序 5中只写出了不同的部分。b由 m 个不同特征向量组成时的 Lanczos收敛情况如图 7所示。图7 m 取不同值时的收敛步数见表 3。表3m1

24、05065803008001000迭代步数6868129169188192可见，随着 m 的增大，迭代步数逐渐增加。理论上，当 b 仅由 A 的 m 个不同特征向量 的线性组合表示时， Lanczos 方法必至多 m 步收敛。但实际上，当 m 较大时，由于精度等 方面的限制， m 个特征向量并不都线性无关，所以，当 m 较大时，迭代步数可能比 m 小的 多，比如本例中 m>300 后，迭代步数远小于 m。而同样由于计算精度限制，在 m 较小时， 迭代步数可能溢出 m 步，比如本例 m=65 和 80 的情形。另外， m 值大于一定值后，随着 m 的增加，收敛步数逐渐趋于定值，可能的原因是

25、： 随着 m 的增大，特征向量线性相关性逐渐增强。附 matlab 程序 4 ： clear all 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);for i=1:7M=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 解法 = ticy=norm(r0); F(1)=no

26、rm(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);%给各变量赋初始值for j=2:n r0=A*Q(:,j)-bita(j-1)*Q(:,j-1);alpha(j)=Q(:,j)'*r0; r0=r0-alpha(j)*Q(:,j);if norm(r0)>1bita(j)=norm(r0); Q(:,j+1)=r0/bita(j);endfor k=1:jT(k,k)=alpha(k);

27、endfor k=1:j-1T(k+1,k)=bita(k);T(k,k+1)=bita(k);%生成三对角阵 Tende(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 及 rF(j,i)=r;if r<sqrt(eps) break ;end %r 的精度达到要求后停止迭代，得到最终 X end toc jclear e y1 X r Tendsemilogy(F(:,1),'r' );hold on ;semilogy(F(:,2),'b'

28、);hold on ;semilogy(F(:,3),'y' );hold on ;semilogy(F(:,4),'black' );hold on;semilogy(F(:,5),'g' );hold on ;semilogy(F(:,6),'m' );hold on ;semilogy(F(:,7),'c' );hold on ;legend( 'm=10' ,'m=50' , 'm=65' , 'm=80' , 'm=300'

29、, 'm=800'title( 'b 由不同特征向量个数组成时 Lanczos 的收敛性 ' )xlabel( ' 迭代步数')ylabel( ' 残差对数log(|rk|)' ), 'm=1000' )第五题构造对称不定矩阵， 验证 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通过以上实验可以看到， 当负特