整理高等数值分析作业第二次实验.docx
《整理高等数值分析作业第二次实验.docx》由会员分享，可在线阅读，更多相关《整理高等数值分析作业第二次实验.docx（23页珍藏版）》请在冰豆网上搜索。
整理高等数值分析作业第二次实验
（二）环境影响经济损益分析的步骤
（3）评价单元划分应考虑安全预评价的特点，以自然条件、基本工艺条件、危险、有害因素分布及状况便于实施评价为原则进行。
2）规划实施可能对环境和人群健康产生的长远影响。
专项规划中的指导性规划 环境影响篇章或说明
2.量化环境影响后果
为了有别于传统的忽视环境价值的理论和方法，环境经济学家把环境的价值称为总经济价值（TEV），包括环境的使用价值和非使用价值两个部分。
（五）安全预评价方法
（6）环境影响评价结论的科学性。
环境影响经济损益分析一般按以下四个步骤进行：
内涵资产定价法基于这样一种理论，即人们赋予环境的价值可以从他们购买的具有环境属性的商品的价格中推断出来。
高等数值分析第二次实验作业
注：
矩阵阶数均为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&jj=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&jj=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&jj=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&jj=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;
whileje
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;
whileje
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
升级会员 