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

《高等数值分析作业-第二次实验》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、高等数值分析第二次实验作业注：矩阵阶数均为1000阶T1.构造例子特征值全部在右半平面时，观察基本的Arnoldi方法和GMRES方法的数值性态，和相应重新启动算法的收敛性。Answer：Ø关于特征值均在右半平面的矩阵A：首先构造对角矩阵D，其中D矩阵的对角是由以下形式的2×2的子矩阵组成其对角元：S=a-bba，其中a，b为实数此时，S矩阵的特征值分别为a+bi和a-bi，这样D=diag(S1,S2,S3……Sn)矩阵的特征值均分布在右半平面。另矩阵A=QTAQ，则A矩阵的特征值也均在右半平面，生成矩阵的过程代码如下所示：N=500%生成的矩阵为2N阶A=ze

2、ros(2*N);%delta控制特征值分布的密集程度，即控制条件数大小，delta越大条件数越大delta=0.1;%构造D矩阵forj=1:NA(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;endU=orth(rand(2*N,2*N));A=U'*A*U;Ø首先进行观察基本的Arnoldi方法和GMRES方法的数值性态，考虑N=500，即矩阵为1000阶，取X0=zeros(N,1)，b=ones(N,1)，收敛准则为e=

3、10-6;Arnoldi方法函数如下：第

4、页13function[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;whilenorm(r)>e&j

5、=1;be1=belta*e1;try第

6、页13[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;enderror(j)=log10(norm(r));endend第

7、页13GMRES方法的函数如下：第

8、页13function[xm,error,num]=GMRES(A,x0,b,e)%ARNOLDISummaryofthisfunctiongoeshere%DetailedexplanationgoeshereN=size(A,1);r

9、=b-A*x0;belta=norm(r);v=r/belta;V=v;j=0;er=1000;whileer>e&j

10、j+1);enderror(j)=log10(er);endend第

11、页13其中UX，LX和minresYK是自己编写的求解（上、下）三角矩阵的函数，代码如下：第

12、页13functionyk=LX(TK,b)n=size(TK,1);yk=zeros(n,1);yk(1)=b(1)/TK(1,1);fori=2:nyk(i)=b(i);forj=1:i-1yk(i)=yk(i)-TK(i,j)*yk(j);endyk(i)=yk(i)/TK(i,i);endendfunctionyk=UX(TK,b)n=size(TK,2);yk=zeros(n,1);yk(n)

13、=b(n)/TK(n,n);fori=1:n-1k=n-i;yk(k)=b(k);forj=k+1:nyk(k)=yk(k)-TK(k,j)*yk(j);endyk(k)=yk(k)/TK(k,k);endendfunctionyk=minresYK(TK,b)n=size(TK,2);yk=zeros(n,1);yk(n)=b(n)/TK(n,n);fori=1:n-1k=n-i;yk(k)=b(k);forj=k+1:nyk(k)=yk(k)-TK(k,j)*yk(j);endyk(k)=yk(k)/TK(k,k);endend第

14、页13两种方法的计算结果如下

15、所示：第

16、