整理高等数值分析作业第二次实验Word文件下载.docx
《整理高等数值分析作业第二次实验Word文件下载.docx》由会员分享,可在线阅读,更多相关《整理高等数值分析作业第二次实验Word文件下载.docx(23页珍藏版)》请在冰豆网上搜索。
关于特征值均在右半平面的矩阵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;
error(j)=log10(norm(r));
GMRES方法的函数如下:
function[xm,error,num]=GMRES(A,x0,b,e)
%ARNOLDISummaryofthisfunctiongoeshere
%Detailedexplanationgoeshere
er=1000;
whileer>
e&
[Q,R]=qr(h(1:
j+1,1:
gk=Q'
*belta*eye(j+1,1);
ym=minresYK(R,gk);
er=gk(j+1);
error(j)=log10(er);
其中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);
yk(i)=yk(i)/TK(i,i);
functionyk=UX(TK,b)
n=size(TK,2);
yk(n)=b(n)/TK(n,n);
fori=1:
n-1
k=n-i;
yk(k)=b(k);
forj=k+1:
yk(k)=yk(k)-TK(k,j)*yk(j);
yk(k)=yk(k)/TK(k,k);
functionyk=minresYK(TK,b)
两种方法的计算结果如下所示:
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)
whileabs(norm(r))>
m
error=abs(norm(r));
function[xm,error]=GMRES_m(A,x0,b,e,m)
whileabs(er)>
error=abs(er);
B、调用上面更改后的代码,编写重启的算法代码如下:
function[xm,err,num]=Arnoldim(A,x0,b,e,m)
jmax=N/m;
xm=x0;
error=1000;
whilej<
jmax&
error>
e
[xm,error]=Arnoldi_m(A,xm,b,e,m);
err(j)=log10(error);
log10(error)
function[xm,err,num]=GMRESm(A,x0,b,e,m)
[xm,error]=GMRES_m(A,xm,b,e,m);
计算结果如下所示(计算时取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算法相应地会怎样?
对A的特征值进行平移,由于只要变换A矩阵的实部,这里将上面的矩阵生成做如下改造:
A