北航数值分析大作业程序1.docx
《北航数值分析大作业程序1.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业程序1.docx（14页珍藏版）》请在冰豆网上搜索。
北航数值分析大作业程序1
《数值分析B》
第一次数值分析大作业
院系：
04能源与动力工程学院
姓名：
王开逍
学号：
SY1104207
一．算法设计方案：
1.首先应该创建原矩阵，题目所给的矩阵是大型的带状矩阵，为了节省存储量，矩阵的带外元素不给存储，仅存矩阵的带内元素，即将其进行压缩存储。
2．由于λ1‹λ2‹…‹λ501，所以在以所有特征值建立的数轴上，λ1、λ501位于数轴的两端，两者之一必为按模最大。
利用幂法，可以求出来按模最大的特征值，但该值可能为λ1也可能为λ501；
3．上步求出按模最大的特征值λM后，将原矩阵平移-λM,再利用幂法求一次平移后矩阵的按模最大的特征值，即是数轴上另一端点值λM′。
4．比较λM与λM+λM′的大小，大的为矩阵A的最大特征值，小的为A矩阵的最小的特征值；
5．利用反幂法，求矩阵A的按模最小的特征值。
但是反幂法中要用到线性方程组的求解，而原矩阵A又是带状矩阵，采用LU分解。
所以在这之前要定义一个LU分解子程序，将A矩阵分解为单位下三角矩阵L和上三角矩阵U的乘积。
6．先利用循环求出k从1到39变化的uk的值。
对于每一个uk值，都将压缩后的原始矩阵A的第三行减去相应的uk，形成新的矩阵A（k），然后调用LU分解的子程序，利用反幂法求出与uk最接近的特征值，该特征值等于利用反幂法求出的值与uk的和。
7．A的谱范数条件数等于按模最大的特征值的绝对值与按模最小的特征值的绝对值之比，按模最大的特征值与按模最小的特征值已分别在前面求出。
6．A的行列式的值就是矩阵A进行LU分解后U的对角线元素的乘积。
二．源程序：
该计算程序使用FORTRAN语言编写
modulelinear1!
生成模块，便于主程序中调用
implicitnone
contains
subroutinelu（matrix,us,b,x）!
这个子程序是用来对矩阵进行LU分解的，并且可以求解线性方程AX=B
real（kind=8）matrix（:
:
）,b（:
）,us（:
:
）,x（:
）
integern
real（kind=8）,allocatable:
:
l（:
:
）,y（:
）
real（kind=8）a,c,e,f
integeri,j,k,t
n=size（matrix,1）
allocate（l（n,n）,y（n））
a=0.0
c=0.0
doj=1,n
us（1,j）=matrix（1,j）
enddo
doi=2,n
l（i,1）=matrix（i,1）/us（1,1）
enddo
dok=2,n-1
doj=k,n
doi=k+1,n
dot=1,k-1
a=a+l（k,t）*us（t,j）
c=c+l（i,t）*us（t,k）
enddo
us（k,j）=matrix（k,j）-a
l（i,k）=（matrix（i,k）-c）/us（k,k）!
矩阵A的LU分解结束
a=0.0
c=0.0
enddo
enddo
enddo
a=0.0
dot=1,n-1
a=a+l（n,t）*us（t,n）
enddo
us（n,n）=matrix（n,n）-a
y
（1）=b
（1）
e=0.0
f=0.0
doi=2,n
dot=1,i-1
e=e+l（i,t）*y（t）
enddo
y（i）=b（i）-e
e=0.0
enddo
x（n）=y（n）/us（n,n）
doi=n-1,1,-1
dot=i+1,n
f=f+us（i,t）*x（t）
enddo
x（i）=（y（i）-f）/us（i,i）!
求得解向量X（i）
f=0.0
enddo
endsubroutinelu
realfunctionuu（aa,bb）!
这个子程序是用来做两个向量啊a,b的乘积的
real（kind=8）aa（:
）,bb（:
）
integern,i
n=size（aa）
uu=0.0
doi=1,n
uu=uu+aa（i）*bb（i）
enddo
endfunction
subroutineau（matrix,u,c）!
这个子程序是用来计算矩阵和向量相乘的，最后返回一个向量值
real（kind=8）matrix（:
:
）,u（:
）,c（:
）
integern,i,k
n=size（matrix,1）
doi=1,n
c（i）=0.0
enddo
dok=1,n
doi=1,n
c（k）=c（k）+matrix（k,i）*u（i）
enddo
enddo
endsubroutineau
subroutinepositive（matrix,u0,bp）!
这个子程序是用幂法计算矩阵的模最大的特征值
real（kind=8）matrix（:
:
）,u0（:
）,bp
integern,i
real（kind=8）,allocatable:
:
u（:
）,y（:
）,c（:
）
real（kind=8）co,nn
n=size（matrix,1）
allocate（u（n）,y（n）,c（n））
bp=0.0
doi=1,n
u（i）=u0（i）
enddo
dowhile（1>0）
co=bp
nn=uu（u,u）**0.5
doi=1,n
y（i）=u（i）/nn
enddo
callau（matrix,y,c）!
调用矩阵和向量相乘的子程序
doi=1,n
u（i）=c（i）
enddo
bp=uu（y,u）
write（*,"（e20.9）"）bp
if（abs（bp-co）/abs（bp）<=0.1e-11）exit
enddo
endsubroutinepositive
subroutineopposite（matrix,u0,bo）!
这个子程序是用反幂法求出矩阵的模最小的特征值
real（kind=8）matrix（:
:
）,u0（:
）,bo
integern,i
real（kind=8）,allocatable:
:
u（:
）,y（:
）,us（:
:
）,x（:
）
real（kind=8）c,nn
bo=0.0
n=size（matrix,1）
allocate（u（n）,y（n）,us（n,n）,x（n））
doi=1,n
u（i）=u0（i）
enddo
dowhile（1>0）
c=bo
nn=sqrt（uu（u,u））
doi=1,n
y（i）=u（i）/nn
enddo
calllu（matrix,us,y,x）!
调用对矩阵进行LU分解的子程序
doi=1,n
u（i）=x（i）
enddo
bo=uu（y,u）
write（*,"（e20.9）"）bo
if（abs（bo-c）/abs（bo）<=0.1e-3）exit
enddo
endsubroutineopposite
endmodule
programmain!
主程序开始
uselinear1!
调用程序模块linear1
implicitnone
integer,parameter:
:
m=501
integeri,j,k
real（kind=8）a0（m,m）,b0（m,m）,c0（m,m）,us（m,m）
real（kind=8）u0（m）,x（m）,c（m）
real（kind=8）r1,r2,bo,bp
real（kind=8）t1,t501,ts,w,t,cond2,tm
real（kind=8）det
doi=1,m!
对矩阵A进行初始化
doj=1,m
a0（i,j）=0.0
enddo
enddo
doi=1,m
a0（i,i）=（1.64-0.024*i）*sin（0.2*i）-0.64*exp（0.1/i）
enddo
doi=2,m
a0（i,i-1）=0.16
enddo
doi=2,m
a0（i-1,i）=0.16
enddo
doi=3,m
a0（i,i-2）=-0.064
enddo
doi=3,m
a0（i-2,i）=-0.064
enddo
doi=1,m
u0（i）=1.0
enddo
doi=1,m
doj=1,m
b0（i,j）=a0（i,j）
enddo
enddo
callpositive（b0,u0,bp）!
调用幂法子程序求矩阵A的模最大的特征值
tm=bp
r1=bp
write（*,"（'原点平移前矩阵A的模最大的特征值为',e20.9）"）r1
doi=1,m
doj=1,m
c0（i,j）=a0（i,j）
enddo
enddo
doi=1,m
c0（i,i）=a0（i,i）-r1
enddo
callpositive（c0,u0,bp）!
调用幂法子程序求原点平移后矩阵的最大特征值
write（*,"（'原点平移后矩阵（A-uI）的模最大的特征值为',e20.9）"）bp
r2=bp+r1
if（r2>r1）then
t501=r2
t1=r1
else
t501=r1
t1=r2
endif
open（unit=30,file='output/thicknessandf.txt'）!
将程序运行的结果存储在一个TXT的文档里面
write（30,"（'最小的特征值为t1=',e20.9,'最大的特征值为t501=',e20.9）"）t1,t501
write（*,"（'最小的特征值为t1=',e20.9,'最大的特征值为t501=',e20.9）"）t1,t501
callopposite（a0,u0,bo）!
调用反幂法子程序，求矩阵A按模最小的特征值
ts=bo**（0.0-1.0）
write（30,"（'模最小的特征值为',e20.9）"）ts
write（*,"（'模最小的特征值为',e20.9）"）ts
dok=1,39
w=t1+k*（t501-t1）/40.0
doi=1,m
c0（i,i）=a0（i,i）-w
enddo
callopposite（c0,u0,bo）!
对于每一个u（k）调用反幂法，求得原点平移后矩阵模最小特征值
t=bo**（0.0-1.0）+w!
原点平移后矩阵模最小特征值加上平移量就可以得到所要求特征值
write（30,"（'与U',i2,'=',e20.9,'最接近的特征值为',e20.9）"）k,w,t
write（*,"（'与U',i2,'=',e20.9,'最接近的特征值为',e20.9）"）k,w,t
enddo
cond2=abs（tm/ts）!
求矩阵的二范数
write（30,"（'矩阵A的条件数cond（A）2=',e20.9）"）cond2
write（*,"（'矩阵A的条件数cond（A）2=',e20.9）"）cond2
calllu（a0,us,u0,x）
det=1.0
doi=1,m
det=det*us（i,i）!
求矩阵的特征值
enddo
write（30,"（'矩阵A的行列式DET（A）=',e20.9）"）det
write（*,"（'矩阵A的行列式DET（A）=',e20.9）"）det
stop
endprogram!
主程序结束
三．输出结果：
最小的特征值为t1=-0.107000933E+02最大的特征值为t501=0.972455120E+01
模最小的特征值为-0.555854581E-02
与U1=-0.101894772E+02最接近的特征值为-0.101829363E+02
与U2=-0.967886105E+01最接近的特征值为-0.958562413E+01
与U3=-0.916824493E+01最接近的特征值为-0.917267607E+01
与U4=-0.865762882E+01最接近的特征值为-0.865228622E+01
与U5=-0.814701271E+01最接近的特征值为-0.809346541E+01
与U6=-0.763639660E+01最接近的特征值为-0.765940833E+01
与U7=-0.712578049E+01最接近的特征值为-0.711968532E+01
与U8=-0.661516438E+01最接近的特征值为-0.661175517E+01
与U9=-0.610454826E+01最接近的特征值为-0.606611000E+01
与U10=-0.559393215E+01最接近的特征值为-0.558510429E+01
与U11=-0.508331604E+01最接近的特征值为-0.511407297E+01
与U12=-0.457269993E+01最接近的特征值为-0.457887511E+01
与U13=-0.406208382E+01最接近的特征值为-0.409651422E+01
与U14=-0.355146770E+01最接近的特征值为-0.355421889E+01
与U15=-0.304085159E+01最接近的特征值为-0.304108999E+01
与U16=-0.253023548E+01最接近的特征值为-0.253397675E+01
与U17=-0.201961937E+01最接近的特征值为-0.200322594E+01
与U18=-0.150900326E+01最接近的特征值为-0.150355918E+01
与U19=-0.998387146E+00最接近的特征值为-0.993558698E+00
与U20=-0.487771034E+00最接近的特征值为-0.487042661E+00
与U21=0.228450775E-01最接近的特征值为0.223074975E-01
与U22=0.533461189E+00最接近的特征值为0.532417509E+00
与U23=0.104407730E+01最接近的特征值为0.105290380E+01
与U24=0.155469341E+01最接近的特征值为0.158944712E+01
与U25=0.206530952E+01最接近的特征值为0.206034137E+01
与U26=0.257592564E+01最接近的特征值为0.255807431E+01
与U27=0.308654175E+01最接近的特征值为0.308024567E+01
与U28=0.359715786E+01最接近的特征值为0.361362115E+01
与U29=0.410777397E+01最接近的特征值为0.409136867E+01
与U30=0.461839008E+01最接近的特征值为0.460302209E+01
与U31=0.512900620E+01最接近的特征值为0.513292891E+01
与U32=0.563962231E+01最接近的特征值为0.559490416E+01
与U33=0.615023842E+01最接近的特征值为0.608092462E+01
与U34=0.666085453E+01最接近的特征值为0.668036226E+01
与U35=0.717147064E+01最接近的特征值为0.729389173E+01
与U36=0.768208675E+01最接近的特征值为0.771711848E+01
与U37=0.819270287E+01最接近的特征值为0.822521984E+01
与U38=0.870331898E+01最接近的特征值为0.864865720E+01
与U39=0.921393509E+01最接近的特征值为0.925419879E+01
矩阵A的条件数cond（A）2=0.192498067E+04
矩阵A的行列式DET（A）=0.277059095E+119
四．讨论迭代初始向量的选取对于计算结果的影响：
1.影响迭代速度。
因为如果初始向量选取不当，有可能使得迭代开始的时刻有可能不同，因此收敛的速度不一样，达到规定精度所需的迭代次数也就不一样。
当迭代次数过大时，计算机可能不支持从而计算不出结果。
2.影响计算结果的精度。
因为选取不同的迭代初始向量，在规定的迭代次数下，计算机计算出的结果的精度是不一样的，有时可能因为误差太大得到错误的结果。
例如当该程序初始向量选择{1，0，0，……0}时得到的结果和{1,1……1}；得到的结果是不一样的，同一迭代次数下，后者的精度更高。
因此一般情况下，我们在给定初始向量时都选择后者。
程序运行结果显示：
升级会员 