数值分析计算实习题目一.docx
《数值分析计算实习题目一.docx》由会员分享，可在线阅读，更多相关《数值分析计算实习题目一.docx（14页珍藏版）》请在冰点文库上搜索。
数值分析计算实习题目一
《数值分析》大作业
（一）
SY1616666XX
设有
的矩阵
其中
矩阵
的特征值
满足
试求:
1.
的值
2.
的与数
最接近的特征值
3.
的（谱范数）条件数
和行列式
.
1.算法设计方案
本题的核心算法是幂法、带原点平移的幂法、反幂法、带原点平移的反幂法和LU分解法，要点在于选择算法时，应使
的所有零元素都不存储。
故算法设计的思路如下：
第一步，对
使用幂法（Powerlaw），可得
的按模最大的特征值，记为
；如果
则
，对
使用带有原点平移的幂法，令平移量
，可得另一端点的特征值
；否则，
，对
使用带有原点平移的幂法，令平移量
，可得另一端点的特征值
；
第二步，对
使用反幂法（InversePower），可得
的按模最小的特征值
（使用LU的Doolittle分解法）
第三步，根据
计算出
然后利用带有原点平移的反幂法求得
其中平移量
反幂法运算39次,可得
;
第四步,根据定义,非奇异的实对称矩阵
的谱范数条件数
其中
分别是矩阵
的模为最大和模为最小的特征值，对于本题，则有
；
然后，由LU分解可知，
，可得
。
由于题目要求算法中应使
的所有零元素均不存储，故构造一个数组C存放
中不为零的元素。
2.全部的源程序
moduleglobal
integer（4）,parameter:
:
n=501,r=2,s=2,m=5
real（8）,parameter:
:
epsilon=1d-12
real（8）,dimension（m,n）:
:
C,C_0=0
real（8）,dimension（n）:
:
u,y,box
real（8）,dimension（39）:
:
mu_k,lambda_ik
real（8）:
:
lambda_1,lambda_501,lambda_s,cond_A2,det_A
end
!
--------------------------------------主程序---------------------------------------------!
PROGRAMMAIN
useglobal
implicitnone
integer（4）:
:
i,j,k
real（8）:
:
eta,beta,beta_1
C_0（1,3:
501）=-0.064d0
C_0（5,1:
499）=-0.064d0
C_0（2,2:
501）=0.16d0
C_0（4,1:
500）=0.16d0
doi=1,n
C_0（3,i）=（1.64-0.024*i）*sin（0.2*i）-0.64*exp（0.1/i）
Enddo
!
===================求解λ_1和λ_501===================!
C=C_0
callPowerlaw（beta）
if（beta<0）then
lambda_1=beta
write（*,12）,'lambda_1','=',lambda_1
C（3,1:
n）=C（3,1:
n）-（lambda_1-1d0）
callPowerlaw（beta）
lambda_501=beta+（lambda_1-1d0）
write（*,12）,'lambda_501','=',lambda_501
else
lambda_501=beta
write（*,12）,'lambda_501','=',lambda_501
C（3,1:
n）=C（3,1:
n）-（lambda_501+1d0）
callPowerlaw（beta）
lambda_1=beta+（lambda_501+1d0）
write（*,12）,'lambda_1','=',lambda_1
endif
!
=====================求解λ_s====================!
C=C_0
callInversePower（beta）
lambda_s=1/beta
write（*,12）,'lambda_s','=',lambda_s
!
===================求解与μ_k相近的λ_ik==================!
dok=1,39
C=C_0
mu_k（k）=lambda_1+k*（lambda_501-lambda_1）/40
C（3,1:
n）=C（3,1:
n）-mu_k（k）
callInversePower（beta）
lambda_ik（k）=1/beta+mu_k（k）
write（*,10）,'lambda_i',k,'=',lambda_ik（k）
enddo
!
===================求解cond（A）2和detA====================!
cond_A2=max（abs（lambda_1）,abs（lambda_501））/abs（lambda_s）
write（*,12）,'cond_A2','=',cond_A2
C=C_0
callDoolittle（）
det_A=product（C（3,1:
n））
write（*,12）,'det_A','=',det_A
10format（1x,a8,i2,a1,es20.12）
12format（1x,a10,a1,es20.12）
end
!
--------------------------------------主程序----------------------------------------!
!
=================幂法例行子程序==================!
subroutinePowerlaw（beta）
useglobal
implicitnone
integer（4）:
:
i,j,k
real（8）:
:
eta,beta,beta_1
beta_1=100d0
k=0
u
（1）=1d0
u（2:
n）=1d0
dowhile（k<=10000）
k=k+1
eta=（dot_product（u（1:
n）,u（1:
n）））**0.5d0
y（1:
n）=u（1:
n）/eta
u
（1）=dot_product（y（1:
3）,C（3:
5,1））
u
（2）=dot_product（y（1:
4）,C（2:
5,2））
doi=r+1,n-r
u（i）=dot_product（y（i-r:
i+r）,C（1:
m,i））
enddo
u（n-1）=dot_product（y（n-3:
n）,C（1:
4,n-1））
u（n）=dot_product（y（n-2:
n）,C（1:
3,n））
beta=dot_product（y（1:
n）,u（1:
n））
if（abs（（beta-beta_1）/beta）<=epsilon）then
exit
else
beta_1=beta
endif
enddo
end
!
================Doolittle分解例行子程序================!
subroutineDoolittle（）
useglobal
implicitnone
integer（4）:
:
i,j,k,t
real（8）:
:
sum
dok=1,n
doj=k,min（（k+s）,n）
sum=0d0
dot=max（1,（k-r）,（j-s））,（k-1）
sum=sum+C（（k-t+s+1）,t）*C（（t-j+s+1）,j）
enddo
C（（k-j+s+1）,j）=C（（k-j+s+1）,j）-sum
enddo
doi=（k+1）,min（（k+r）,n）
if（k>=n）then
exit
endif
sum=0d0
dot=max（1,（i-r）,（k-s））,（k-1）
sum=sum+C（（i-t+s+1）,t）*C（（t-k+s+1）,k）
enddo
C（（i-k+s+1）,k）=（C（（i-k+s+1）,k）-sum）/C（（s+1）,k）
enddo
enddo
end
!
=================反幂法例行子程序==================!
subroutineInversePower（beta）
useglobal
implicitnone
integer（4）:
:
i,j,k
real（8）:
:
eta,beta,beta_1
beta_1=100d0
k=0
callDoolittle（）
u
（1）=1d0
u（2:
n）=1d0
dowhile（k<=10000）
k=k+1
eta=（dot_product（u（1:
n）,u（1:
n）））**0.5d0
y（1:
n）=u（1:
n）/eta
box=y
callSolve_mu_k（）
beta=dot_product（box（1:
n）,u（1:
n））
if（abs（（1d0/beta-1d0/beta_1）*beta）<=epsilon）then
exit
else
beta_1=beta
endif
enddo
end
!
================求解u_k例行子程序==================!
subroutineSolve_mu_k（）
useglobal
implicitnone
integer（4）:
:
i,j,k,t
real（8）:
:
sum
doi=2,n
sum=0d0
dot=max（1,（i-r））,（i-1）
sum=sum+C（（i-t+s+1）,t）*y（t）
enddo
y（i）=y（i）-sum
enddo
u（n）=y（n）/C（（s+1）,n）
doi=（n-1）,1,-1
sum=0d0
dot=（i+1）,min（n,（i+s））
sum=sum+C（（i-t+s+1）,t）*u（t）
enddo
u（i）=（y（i）-sum）/C（（s+1）,i）
enddo
end
3.计算结果
计算结果如图1所示
图1计算结果
4.讨论迭代初始向量的选取对计算结果的影响,并说明原因。
在MAIN程序中,为初始向量u赋初值
得到上述结果。
下面改变初始向量u的赋值,观察计算结果的差异。
时,结果如图2
图2
时,结果如图3
图3
时,结果如图4
图4
时,结果如图5
图5
时,结果如图6
图6
时,结果如图7
图7
通过观察发现,矩阵迭代法的初始向量选取非常重要,并不是任意选择的初始向量都能收敛得到想要的结果,不同的初始向量可能是计算结果收敛于不同阶的特征值与特征向量,而不一定收敛于第一阶.根据文献1的算例,
其特征值显然为：
相应的特征向量为：
通过不同的初始向量得到的迭代结果如表1。
表1不同初始向量迭代结果
给定的初始
向量的转置
计算得到的特征值
计算得到的特征向量的转置
1（1,1,1,1）
4.000
（1.00000，0.00002，0.00000，0.00000）
2（0,1,1,1）
3.000
（0.00000，1.00000，0.00002，0.00000）
3（0.0001,1,0,0）
4.000
（1.00000，0.00002，0.00000，0.00000）
4（0.00006,1,0,0）
4.000
（1.00000，0.00003，0.00000，0.00000）
5（0.00003,1,0,0）
3.000
（0.00004，1.00000，0.00002，0.00000）
6（0.00004,1,0,0）
3.000
（1.00000，0.00003，0.00000，0.00000）
7（0.00003,1,1,1）
4.000
（1.00000，0.00003，0.00000，0.00000）
从表1中可以看到:
由于选取初始向量的不同,矩阵迭代的结果不一定收敛于X1与λ1,
如表中②、⑤栏所示;若初始向量很接近于某一阶特征向量,如③、④、⑤、⑥栏中的情
况,则收敛结果可能是最接近的这阶向量,也可能收敛于其他阶特征向量。
如给定的向量恰使a1＝0,则迭代的结果将得不到λ1与X1（如表1中第②栏）。
但是,由于事先不知道特征向量,所以无法确定a1是否为零。
实际计算中并不一定需要取n个值,分别取几个线性无关的向量作为初始向量计算找出其中的最大特征值即可,因为一般不会凑巧每一个初始向量都使
。
矩阵迭代法求矩阵的特征值与特征向量时,总可以得到收敛的值,但是收敛于哪一阶却与初始向量的选取很有关系,得到的可能不是第一阶特征值与特征向量。
当选取n个线性无关的向量作为初始向量进行n次计算时,总可以得到它的第一阶特征值与特征向量。
升级会员 