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

高等数值分析第一次实验第一题：构造例子说明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, 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：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=0;semilogy(0,norm(r),r-);hold on;while k1000 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(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 k1000 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;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 k1000 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-); hold on; k=k+1;endtitle(条件数的大小对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：clear allclc%条件数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-);hold on;while k1000 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=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;%初始搜索方向k=0;semilogy(0,norm(r),b-);hold on;while k1000 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(迭代步数)ylabel(残差对数log(|rk|) 第二题对于同样的例子，比较CG和Lanczos的计算结果解：采用和第一题相同的构造方法，构造三个1000阶正定对称矩阵，使条件数k分别为：1000,100,10。分别采用CG和Lanczos方法计算Ax=b，且都设置停机准则为：norm(rk)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)=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=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(迭代步数)ylabel(残差对数log(|rk|) 第三题当A只有m个不同特征值时，对于大的m和小的m，观察有限精度下Lanczos方法如何收敛。解： 分别构建m = 10、50、100、1000四个矩阵A，设置停机准则为：norm(rk)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)=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)300后，迭代步数远小于m。而同样由于计算精度限制，在m较小时，迭代步数可能溢出m步，比如本例m=65和80的情形。另外，m值大于一定值后，随着m的增加，收敛步数逐渐趋于定值，可能的原因是：随着m的增大，特征向量线性相关性逐渐增强。附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); 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)=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 rsqrt(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),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由不同特征向量个数组成时Lanczos的收敛性)xlabel(迭代步数)ylabel(残差对数log(|rk|)第五题构造对称不定矩阵，验证Lanczos方法的近似中断，观察收敛曲线中的峰点个数和特征值的分布关系；观测当出现峰点时，MINRES方法的收敛形态怎样。解：分别构建负特征值个数为m = 0、10、100、500矩阵1000阶系数矩阵A，用Lanczos算法和MINRES算法求解Ax = b，设