高等数值分析作-第二次实验.docx

第 | 页 PAGE \* MERGEFORMAT1 高等数值分析第二次实验作业 注：矩阵阶数均为1000阶 T1. 构造例子特征值全部在右半平面时，观察基本的Arnoldi方法和GMRES方法的数值性态，和相应重新启动算法的收敛性。 Answer： 关于特征值均在右半平面的矩阵A： 首先构造对角矩阵D，其中D矩阵的对角是由以下形式的2×2的子矩阵组成其对角元： S= 此时，S矩阵的特征值分别为a+bi和a-bi，这样D=diag(S1,S2,S3……Sn)矩阵的特征值均分布在右半平面。另矩阵A=QTAQ，则A矩阵的特征值也均在右半平面，生成矩阵的过程代码如下所示： N=500 %生成的矩阵为2N阶 A=zeros(2*N); %delta控制特征值分布的密集程度，即控制条件数大小，delta越大条件数越大 delta=0.1; %构造D矩阵 for j=1:N A(2*j-1,2*j-1)=N+j*delta; A(2*j-1,2*j)=-N-j*delta; A(2*j,2*j-1)=N+j*delta; A(2*j,2*j)=N+j*delta; end U = orth(rand(2*N,2*N)); A = U*A*U; 首先进行观察基本的Arnoldi方法和GMRES方法的数值性态，考虑N=500，即矩阵为1000阶，取X0=zeros(N,1)，b=ones(N,1)，收敛准则为e=10-6; Arnoldi方法函数如下： function [xm,error,num]=Arnoldi(A,x0,b,e) N=size(A,1); r=b-A*x0; belta=norm(r); v=r/belta; V=v; j=0; while norm(r)e jN j=j+1; num(j)=j; w=A*V(1:N,j); for i=1:j h(i,j)=V(1:N,i)*w; w=w-h(i,j)*V(1:N,i); end h(j+1,j)=norm(w); if h(j+1,j) ~=0 v=w/h(j+1,j); V=[V v]; end e1=zeros(j,1); e1(1)=1; be1=belta*e1; try [L,U,S]=lu(h(1:j,1:j)); be1=S*be1; lym=LX(L,be1); ym=UX(U,lym); xm=x0+V(1:N,1:j)*ym; r=b-A*xm; end error(j)=log10(norm(r)); end end GMRES方法的函数如下： function [xm,error,num] = GMRES(A,x0,b,e) %ARNOLDI Summary of this function goes here % Detailed explanation goes here N=size(A,1); r=b-A*x0; belta=norm(r); v=r/belta; V=v; j=0; er=1000; while er e jN j=j+1; num(j)=j; w=A*V(1:N,j); for i=1:j h(i,j)=V(1:N,i)*w; w=w-h(i,j)*V(1:N,i); end h(j+1,j)=norm(w); if h(j+1,j) ~=0 v=w/h(j+1,j); V=[V v]; end try [Q,R]=qr(h(1:j+1,1:j)); gk=Q*belta*eye(j+1,1); ym=minresYK(R,gk); xm=x0+V(1:N,1:j)*ym; er=gk(j+1); end error(j)=log10(er