高等数值分析作业第二次实验.docx
- 文档编号:29195318
- 上传时间:2023-07-21
- 格式:DOCX
- 页数:24
- 大小:418.76KB
高等数值分析作业第二次实验.docx
《高等数值分析作业第二次实验.docx》由会员分享,可在线阅读,更多相关《高等数值分析作业第二次实验.docx(24页珍藏版)》请在冰豆网上搜索。
高等数值分析作业第二次实验
高等数值分析第二次实验作业
注:
矩阵阶数均为1000阶
T1.构造例子特征值全部在右半平面时,观察基本的Arnoldi方法和GMRES方法的数值性态,和相应重新启动算法的收敛性。
Answer:
Ø关于特征值均在右半平面的矩阵A:
首先构造对角矩阵D,其中D矩阵的对角是由以下形式的2×2的子矩阵组成其对角元:
此时,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矩阵
forj=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;
whilenorm(r)>e&j j=j+1; num(j)=j; w=A*V(1: N,j); fori=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); ifh(j+1,j)~=0 v=w/h(j+1,j); V=[Vv]; 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) %ARNOLDISummaryofthisfunctiongoeshere %Detailedexplanationgoeshere N=size(A,1); r=b-A*x0; belta=norm(r); v=r/belta; V=v; j=0; er=1000; whileer>e&j j=j+1; num(j)=j; w=A*V(1: N,j); fori=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); ifh(j+1,j)~=0 v=w/h(j+1,j); V=[Vv]; 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); end end 其中UX,LX和minresYK是自己编写的求解(上、下)三角矩阵的函数,代码如下: functionyk=LX(TK,b) n=size(TK,1); yk=zeros(n,1); yk (1)=b (1)/TK(1,1); fori=2: n yk(i)=b(i); forj=1: i-1 yk(i)=yk(i)-TK(i,j)*yk(j); end yk(i)=yk(i)/TK(i,i); end end functionyk=UX(TK,b) n=size(TK,2); yk=zeros(n,1); yk(n)=b(n)/TK(n,n); fori=1: n-1 k=n-i; yk(k)=b(k); forj=k+1: n yk(k)=yk(k)-TK(k,j)*yk(j); end yk(k)=yk(k)/TK(k,k); end end functionyk=minresYK(TK,b) n=size(TK,2); yk=zeros(n,1); yk(n)=b(n)/TK(n,n); fori=1: n-1 k=n-i; yk(k)=b(k); forj=k+1: n yk(k)=yk(k)-TK(k,j)*yk(j); end yk(k)=yk(k)/TK(k,k); end end 两种方法的计算结果如下所示: Delta=1即条件数较大 Delta=0.1即条件数较小 可以看到delta=0.1时两种方法收敛都很快,delta=1时,两种方法略有差别,GMRES方法收敛速度略快一些,下面就去delta=1的情况进行第一题的讨论。 delta=1时,计算结果如下: 方法 Arnoldi GMRES方法 -6.4922 -6.0110 迭代次数 25 26 运行时间(s) 0.109216 0.135968 可以看到两种方法的迭代次数差不多,有曲线可以看到GMRES方法的的性态较Arnoldi方法良好,Arnoldi方法会有平台段,但是GMRES方法的则平稳地收敛,这也是最后迭代次数GMRES较Arnoldi方法少一次的原因。 Ø对于重启的Arnoldi方法和GMRES方法代码如下所示: A、首先对之前的标准Arnoldi方法和GMRES方法进行修改,代码如下(m为间隔m步重启): function[xm,error]=Arnoldi_m(A,x0,b,e,m) N=size(A,1); r=b-A*x0; belta=norm(r); v=r/belta; V=v; j=0; whileabs(norm(r))>e&j j=j+1; w=A*V(1: N,j); fori=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); ifh(j+1,j)~=0 v=w/h(j+1,j); V=[Vv]; 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 end error=abs(norm(r)); end function[xm,error]=GMRES_m(A,x0,b,e,m) N=size(A,1); r=b-A*x0; belta=norm(r); v=r/belta; V=v; j=0; er=1000; whileabs(er)>e&j j=j+1; w=A*V(1: N,j); fori=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); ifh(j+1,j)~=0 v=w/h(j+1,j); V=[Vv]; 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=abs(er); end end B、调用上面更改后的代码,编写重启的算法代码如下: function[xm,err,num]=Arnoldim(A,x0,b,e,m) N=size(A,1); jmax=N/m; j=0; xm=x0; error=1000; whilej j=j+1; [xm,error]=Arnoldi_m(A,xm,b,e,m); num(j)=j; err(j)=log10(error); end j log10(error) end function[xm,err,num]=GMRESm(A,x0,b,e,m) N=size(A,1); jmax=N/m; j=0; xm=x0; error=1000; whilej j=j+1; [xm,error]=GMRES_m(A,xm,b,e,m); num(j)=j; err(j)=log10(error); end j log10(error) end 计算结果如下所示(计算时取m=5,即每隔5步进行算法的重启): delta=1m=5(每隔5步重启) delta=1m=10(每隔10步重启) delta=10m=5(每隔5步重启) delta=10m=10(每隔10步重启) 当delta=100m=5(每隔5步重启),结果如下: 从结果可以看出: 1.重启算法的总迭代次数(重启次数×m)要多于基本的算法; 2.重启的Arnoldi方法收敛速度要比GMRES方法快; 3.对于这个问题,重启算法也可以稳定地收敛,及每次重启后的x0都能有所改善; 4.当条件数增大时,明显重启次数会增加; 5.总的来说,对题中的矩阵A,四种方法均能稳定地收敛。 T2.对于1中的矩阵,将特征值进行平移,使得实部有正有负,和1的结果进行比较,方法的收敛速度会如何? 基本的Arnoldi算法有无峰点? 若有,基本的GMRES算法相应地会怎样? Answer: Ø对A的特征值进行平移,由于只要变换A矩阵的实部,这里将上面的矩阵生成做如下改造: A=zeros(2*N); delta=0.1; s=1.5%s用于控制实部为负的特征值的个数,s越大,实部为负的特征值个数越多。 forj=1: N A(2*j-1,2*j-1)=-s*delta+j*delta; A(2*j-1,2*j)=s*delta-j*delta; A(2*j,2*j-1)=-s*delta+j*delta; A(2*j,2*j)=-s*delta+j*delta; end U=orth(rand(2*N,2*N)); A=U'*A*U; Ø基本算法的结果如下所示: 2个实部为负的特征值 10个实部为负的特征值 100个实部为负的特征值 200个实部为负的特征值 500个实部为负的特征值 900个实部为负的特征值 990个实部为负的特征值 998个实部为负的特征值 1000个实部为负的特征值 实部为负的特征值个数 2 10 100 200 500 Arnoldi方法 -6.0416 -6.0025 -11.8957 -12.2160 -12.2839 迭代次数 734 883 1000 1000 1000 运行时间(s) 17.517751 31.238844 41.911802 41.605778 41.496951 GMRES方法 -6.0008 -6.0702 -12.4162 -12.6358 -12.9732 迭代次数 726 880 1000 1000 1000 运行时间(s) 15.142969 27.069932 37.133932 36.937593 36.433673 实部为负的特征值个数 900 990 998 1000 Arnoldi方法 -12.0507 -6.0954 -6.0192 -6.0251 迭代次数 1000 882 738 656 运行时间(s) 42.070828 28.833792 17.203467 12.627750 GMRES方法 -12.5621 -6.0702 -6.0105 -6.0043 迭代次数 1000 878 730 646 运行时间(s) 36.308397 24.999022 16.011876 10.850091 Ø重启算法的计算结果如下所示: 2个实部为负的特征值,每10步重启 2个实部为负的特征值,每100步重启 10个实部为负的特征值,每10步重启 10个实部为负的特征值,每100步重启 100个实部为负的特征值,每10步重启 100个实部为负的特征值,每100步重启 500个实部为负的特征值,每10步重启 500个实部为负的特征值,每100步重启 998个实部为负的特征值,每10步重启 500个实部为负的特征值,每100步重启 Ø结论: 1.对于上面的矩阵,基本的Arnoldi方法与GMRES方法均能收敛,通知GMRES方法收敛速度相对较快,而重启的算法则均不收敛; 2.可以看到单特征值均匀地分布在左右两边平面时(即有一半的特征值实部为负,另一半实部为正),问题的病态程度最高; 3.问题越病态,重启的Arnoldi算法得到的结果越不好,而重启的GMRES对于结果基本就没有改进,这两种方法均不收敛; 4.基本的Arnoldi算法有峰点,这是由于实部为负的特征值的存在,使得H矩阵发生近似奇异导致发生近似中断而引起的,而GMRES算法则相对稳定,当然这两种算法都能收敛。 5.Arnoldi方法的峰点个数与其特征值的分布有关系,特征在左右半平面分布越均匀,峰点个数越多,当实部为负数的特征值小于N/2(N为矩阵阶数),峰点个数等于这类特征值个数,当实部为负数的特征值大于N/2(N为矩阵阶数),峰点个数等于N减这类特征值个数。 6.两种基本算法收敛速度在前期均很慢,但是到迭代后期收敛速度突增,同时可以看到,当问题相当病态时(例如上面100个、200个、500个实部为负的特征值对应的情况),两种方法均是最后一步才一下子收敛。 T3.对1中的例子固定特征值的实部,变化虚部,比较收敛性。 Answer: Ø首先对T1的代码进行简单修改,实部不变,虚部进行一定的放大(引入系数α),相应代码如下: N=500 A=zeros(2*N); x0=zeros(2*N,1); b=ones(2*N,1); e=10^(-6); delta=1; alpha=1; forj=1: N A(2*j-1,2*j-1)=N+j*delta; A(2*j-1,2*j)=-N-alpha*j*delta; A(2*j,2*j-1)=N+alpha*j*delta; A(2*j,2*j)=N+j*delta; end U=orth(rand(2*N,2*N)); A=U'*A*U; 计算结果如下所示,取delta=1,α分别为1、3、5,取X0=zeros(N,1),b=ones(N,1),收敛准则为e=10-6: Arnoldi方法 GMRES 重启Arnoldi方法,每5步重启 重启GMRES方法,每5步重启 重启Arnoldi方法,每10步重启 重启GMRES方法,每10步重启 重启GMRES方法,每10步重启: Ø结论: 1.随着虚部被放大,四种算法的收敛速度均越来越慢; 2.虽然速度变慢但是,两种基本的算法均能收敛, 与迭代次数基本成线性关系,而重启GMRES算法也都能收敛,可以看到即使是α为100的情况,重启GMRES算法每一次重启得到的 随着重启次数基本线性减小,即重启后的x0会有所改善; 3.对于重启Arnoldi算法,他是不稳定的,当α过大,而重启间隔的步数较小时,可能发生不收敛的情况,例如上面的α=5,m-5,重启后的结果可能越来越差,并没有得到改善; T4.当A只有m个不同特征值时,对于大的m和小的m,观察Arnoldi方法和GMRES方法的收敛性. Answer: A的生成参考第一次实验作业的代码,如下所示: N=1000 m=10; DIA=linspace(1,1000,m); VEC=zeros(N,1); k=N/m; i=1; ii=0; whilei<=m forj=1: k VEC(ii+j)=DIA(i); end ii=ii+k; i=i+1; end D=diag(VEC); U=orth(rand(N,N)); A=U'*D*U; 对构造好的A进行计算,在计算时,取X0=zeros(N,1);b=ones(N,1);e=10-6,计算结果如下: Arnoldi方法 GMRES方法 相同的特征值个数 10 20 50 100 500 1000 Arnoldi方法 -10.6206 -10.9011 -6.0109 -6.2994 -6.0594 -6.0024 迭代次数 10 20 42 62 128 165 运行时间(s) 0.021357 0.038060 0.082905 0.132112 0.370289 0.584533 GMRES方法 -10.6400 -11.0604 -6.0084 -6.0891 -6.0663 -6.0048 迭代次数 10 20 42 61 125 162 运行时间(s) 0.013224 0.023461 0.055279 0.093193 0.270723 0.447675 显然,上面构造出来的矩阵的条件数是相同的,得到以下几个结论: 1.对于基本的Arnoldi方法和GMRES方法,计算结果总是收敛的,GMRES方法收敛性较Arnoldi方法好,当m较大时,GMRES总能比Arnoldi方法提前几步收敛; 2.如果A只有m个不同的特征值,两种方法至多m步就可以找到精确解; 3.在m较大的时候,算法收敛所需要的迭代次数远小于m,当m较小时,可能需要接近于m步才能找到准确解; 4.由上面计算可以看到在m值较小时,Arnoldi方法计算出来的第k步的误差 会有一个在计算时会有先增大后减小的过程,最后收敛,而GMRES方法则没有这一性质,结果单调地收敛到精确解; 5.上面的性质与实验1中的Lanczos方法得到的结果类似。 T5.取初始值近似解为零向量,右端项b仅由A的m个不同特征向量的线性组合表示时,Arnoldi方法和GMRES方法的收敛性如何? Answer: Ø同上面那题一样,参考第一次实验的矩阵A生成的方法,我们可以知道U的每一个行向量均为对应A的特征值(D的某一个元素)的特征向量,因此构造m不同时的右端项b代码如下所示: N=1000 m=10; DIA=linspace(1,N,N); D=diag(DIA); U=orth(rand(N,N)); A1=U'*D*U; b1=zeros(N,1); fori=1: m b1=b1+rand (1)*(U(i,: ))'; end X0=zeros(N,1); e=10^(-6); tic [S1,ek,num]=LanczosT3(A1,X0,b1,e); toc 对构造好的A进行计算,在计算时,取X0=zeros(N,1);e=10-6,计算结果如下: Arnoldi方法 GMRES方法 特征向量个数 10 20 50 100 500 750 1000 Arnoldi方法 -6.0597 -6.1419 -6.0282 -6.0041 -6.0778 -6.0089 -6.0549 迭代次数 44 87 118 132 153 155 164 运行时间 0.08051 0.17837 0.290995 0.34232 0.43846 0.45063 0.507530 GMRES 方法 -6.0244 -6.0286 -6.0567 -6.0291 -6.0311 -6.0115 -6.0510 迭代次数 43 88 112 131 148 152 160 运行时间 0.04774 0.13337 0.19433 0.25138 0.32558 0.352942 0.385624 得到以下几个结论: 1.对于基本的Arnoldi方法和GMRES方法,计算结果总是收敛的,GMRES方法收敛性较Arnoldi方法好,当m较大时,GMRES总能比Arnoldi方法提前几步收敛; 2.M较小时,计算过程中收敛误差会产生波动,m较大时,波动减小,且Arnoldi方法波动较GMRES方法剧烈; 3.收敛步数随着m值的增大而增加,但是由上表我们可以看出收敛步数的增幅逐渐减小,可以推测当m达到一定值后,收敛步数可能趋近于某一定值; 4.这个与第一次实验的Lanczos方法的结论类似。 5. 6. 7. 8.(注: 素材和资料部分来自网络,供参考。 请预览后才下载,期待你的好评与关注! ) 9. 10. 11.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 高等 数值 分析 作业 第二次 实验