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

1、高等数值分析第一次实验第一题：构造例子说明CG的数值形态。当步数 = 阶数时CG的解如何？当A的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时，方法的收敛性如何？解：用Housholder变换和对角阵构造1000阶正定对称矩阵A：1） 构造对角阵D = diag( linspace(1, 1000, 1000) )；2） 构造Householder阵H。取单位向量w=1,0,0,.0T，I为1000阶单位矩阵，H = I wTw。3） 构造对称正定矩阵A。A = HTDH。由于D是对角阵，H是对称的，所以A对称；且A与D具有相同的特征值linspace(1, 1000,

2、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)

3、)时，条件数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'A=H'*D*H;%生成1000阶的对称矩阵b=rand(1000,1);x=zeros(1000,1);%初始解r=b-A*x;%初始残量p=r;%初始搜索方向k=

4、0;semilogy(0,norm(r),'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),'r-'); hold on; k=k+1;end%条件数k=100clear allD=diag(linspace(1,100,1000);w=eye(1,1000);I=eye

5、(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<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,

6、norm(r),'b-'); hold on; k=k+1;end%条件数k=10clear 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;%初始残量p=r;%初始搜索方向k=0;semilogy(0,norm(r),'black-');hold on;while k<1000 alpha = r'*p/(p'*A*

7、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-'); hold on; k=k+1;endtitle('条件数的大小对CG法收敛特性的影响');xlabel('迭代步数')ylabel('残差对数log(|rk|)') 2.构造特殊特征值分布构造对称正定矩阵A1和A2。D1=diag( linspace(1, 1000, 1000)

8、)时，条件数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：clear allclc%条件数k=1000,特

9、征值均匀分布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-');hold on;while k<1000 alpha = r'*p/(p'*A*p); x = x+alpha*p; rold = r; r = rold-alpha*A*p; beta

10、= r'*r/(rold'*rold); p = r+beta*p; semilogy(k,norm(r),'r-'); hold on; k=k+1;end%条件数k=1000，最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值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

11、;%初始搜索方向k=0;semilogy(0,norm(r),'b-');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),'b-'); hold on; k=k+1;endtitle('特征值分布对CG法收敛特性的影响');xlabel('迭代步数'

12、)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对比LanczosCG条件数迭代步数时间/s迭代步数时间/s10002430.822670.411001240.431340.2110410.24430

13、.07 可以看到，Lanczos的迭代步数比CG小，说明Lanczos收敛更快，但是其运算时间却明显大于CG，因此总体上CG的收敛效率要比Lanczos高。因此在计算对称正定问题时，优先选用CG算法。附matlab程序2：clear allclcn=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解法 =ti

14、cy=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); %给各变量赋初始值 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)>1 bita(j)=norm(r0); Q(:,j+1)=r0/bita(j);

15、 end for k=1:j T(k,k)=alpha(k); end for k=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; if r/norm(b)<1e-8 break; end %r的精度达到要求后停止迭代，得到最终X end semilogy(F,'r');hold on;j %迭代步数toc%=CG算法=ticp=

16、R0;for i=1:n R=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 break; end beta=(R'*R)/(R0'*R0); p=R+beta*p; R0=R; endsemilogy(G,'b');toci %迭代步数legend('Lanczos','CG')title('Lanczos与CG算法收敛性对比,条件数k=10')xlabel('迭代步数

17、')ylabel('残差对数log(|rk|)') 第三题当A只有m个不同特征值时，对于大的m和小的m，观察有限精度下Lanczos方法如何收敛。解： 分别构建m = 10、50、100、1000四个矩阵A，设置停机准则为：norm(rk)<sqrt(eps)，分别计算Ax = b的解，如图6所示。图 6各个m值达到收敛时迭代步数如表2所示。表2m5008009509901000迭代步数162856119242可见，相同的特征值个数越多，Lanczos收敛速度越快，而且只要m值稍小于矩阵阶数，收敛速度就明显提高，如图中m=990比m=1000的收敛速度快了一倍。事

18、实上，由于所构造的矩阵条件数都不坏（最大为1000），因此算法的收敛步数远小于m值。附matlab程序3：clear allclcfor i=1:5n=1000; M=500,800,950,990,1000;m=M(i);%m = 500、800、950、990、1000D=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;

19、%=Lanczos解法 =ticy=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); %给各变量赋初始值 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)>1 bita(j)=norm(r0); Q(:,j

20、+1)=r0/bita(j); end for k=1:j T(k,k)=alpha(k); end for k=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; if r/norm(b)<sqrt(eps) break; end %r的精度达到要求后停止迭代，得到最终X end tocjclear e y1 W X r Tendsemilog

21、y(F(:,1),'r');hold on;semilogy(F(:,2),'b');hold on;semilogy(F(:,3),'y');hold on;semilogy(F(:,4),'black');hold on;semilogy(F(:,5),'g');hold on;legend('m=500','m=800','m=950','m=990','m=1000')title('不同特征值个数的Lanczos收敛性&

22、#39;)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所示。图 7m取

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

24、附matlab程序4：clear allclcn=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)=norm(r0); Q(:,1)=r0/norm(r0); r0=A*Q(:,1);

25、 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)>1 bita(j)=norm(r0); Q(:,j+1)=r0/bita(j); end for k=1:j T(k,k)=alpha(k); end for k=1:j-1 T(k+1,k)=bi

26、ta(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; if r<sqrt(eps) break; end %r的精度达到要求后停止迭代，得到最终X end tocjclear e y1 X r Tendsemilogy(F(:,1),'r');hold on;semilogy(F(:,2),'b');hold on;semilogy(F(:,3),&

27、#39;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','m=800','m=1000')title('b

28、由不同特征向量个数组成时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通过以上实验可以看到，当负特征值个数不超过特征值总个数一半的时候，负特征值越多，峰点个数越多。实际上，更一般的结论