北航数值分析大作业2Word文档格式.docx
《北航数值分析大作业2Word文档格式.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业2Word文档格式.docx（30页珍藏版）》请在冰豆网上搜索。
这个子程序是用来做两个向量啊a,b的乘积的
real（kind=8）aa（:
）,bb（:
integern,i
n=size（aa）
uu=0.0
uu=uu+aa（i）*bb（i）
endfunction
subroutineuptr（ma,ma1）!
这个子程序是用来对矩阵A进行拟上三角化的
real（kind=8）ma（:
）,ma1（:
integern
real（kind=8）,allocatable:
:
u（:
）,p（:
）,q（:
）,w（:
）,wu（:
）,up（:
）,at（:
real（kind=8）t,d,c0,h,he
integeri,j,r
n=size（ma,1）
allocate（u（n）,p（n）,q（n）,w（n）,wu（n,n）,up（n,n）,at（n,n）,c（n））
doj=1,n
ma1（i,j）=0.0
enddo!
ma1的零初始化
dor=1,n-2
he=0.0
doi=r+1,n
he=he+ma（i,r）**2
d=sqrt（he）
if（ma（r+1,r）>
0）then
c0=0.0-d
else
c0=d
endif
doi=1,r
ma1（i,r）=ma（i,r）
ma1（r+1,r）=c0
doi=r+2,n
if（he==0）then
cycle
h=c0**2-c0*ma（r+1,r）
u（i）=0.0
u（r+1）=ma（r+1,r）-c0
u（i）=ma（i,r）
at（i,j）=ma（j,i）
callau（at,u,c,n）
p（i）=c（i）/h
callau（ma,u,c,n）
q（i）=c（i）/h
t=uu（p,u）/h
w（i）=q（i）-t*u（i）
ma（i,j）=ma（i,j）-w（i）*u（j）-u（i）*p（j）
ma1（i,n-1）=ma（i,n-1）
ma1（i,n）=ma（i,n）
endsubroutineuptr
subroutinetwo（d,s1,s2,a1,a2）!
这个子程序是用来求二阶子阵的两个特征值的
real（kind=8）d（:
complex（kind=8）s1,s2
real（kind=8）a1,a2,b1,b2
n=size（d,1）
if（（d（n-1,n-1）+d（n,n））**2-4*（d（n-1,n-1）*d（n,n）-d（n-1,n）*d（n,n-1））<
0.0）then
a1=（d（n-1,n-1）+d（n,n））/2
a2=（d（n-1,n-1）+d（n,n））/2
b1=sqrt（0.0-（（d（n-1,n-1）+d（n,n））**2-4*（d（n-1,n-1）*d（n,n）-d（n-1,n）*d（n,n-1））））/2
b2=0.0-sqrt（0.0-（（d（n-1,n-1）+d（n,n））**2-4*（d（n-1,n-1）*d（n,n）-d（n-1,n）*d（n,n-1））））/2
a1=（d（n-1,n-1）+d（n,n））/2+sqrt（（d（n-1,n-1）+d（n,n））**2-4*（d（n-1,n-1）*d（n,n）-d（n-1,n）*d（n,n-1）））/2
a2=（d（n-1,n-1）+d（n,n））/2-sqrt（（d（n-1,n-1）+d（n,n））**2-4*（d（n-1,n-1）*d（n,n）-d（n-1,n）*d（n,n-1）））/2
b1=0.0
b2=0.0
s1=cmplx（a1,b1）
s2=cmplx（a2,b2）
endsubroutinetwo
subroutineak2（ak,mk,ak1,m）!
这个子程序是对A（K）进行带双步位移QR分解的方法，并且求出A（K+1）
real（kind=8）ak（:
）,mk（:
）,ak1（:
integerm
real（kind=8）s,t0,c
integeri,j,t
s=ak（m-1,m-1）+ak（m,m）
t0=ak（m-1,m-1）*ak（m,m）-ak（m,m-1）*ak（m-1,m）
doi=1,m
doj=1,m
c=0.0
dot=1,m
c=c+ak（i,t）*ak（t,j）
if（i/=j）then
mk（i,j）=c-s*ak（i,j）
mk（i,j）=c-s*ak（i,j）+t0
callmkak（ak,mk,ak1,m）
endsubroutineak2
subroutinemkak（ak,mk,ak1,m）!
!
用来对中间矩阵M（K）进行QR分解的，并且求出A（K+1）
real（kind=8）mk（:
）,ak（:
real（kind=8）d,c0,h,t,he
b（:
）,c1（:
）,v（:
）,bt（:
）,c1t（:
allocate（b（m,m）,c1（m,m）,u（m）,v（m）,bt（m,m）,c1t（m,m）,p（m）,q（m）,w（m）,c（m））
b（i,j）=mk（i,j）
c1（i,j）=ak（i,j）
dor=1,m-1
doi=r+1,m
he=he+b（i,r）**2
doi=r,m
if（b（r,r）>
h=c0**2-c0*b（r,r）
u（r）=b（r,r）-c0
u（i）=b（i,r）
bt（i,j）=b（j,i）
c1t（i,j）=c1（j,i）
callau（bt,u,c,m）
v（i）=c（i）/h
b（i,j）=b（i,j）-u（i）*v（j）
callau（c1t,u,c,m）
callau（c1,u,c,m）
c1（i,j）=c1（i,j）-w（i）*u（j）-u（i）*p（j）
ak1（i,j）=c1（i,j）
endsubroutinemkak
subroutineqr（ma,ma1,q1,rq）!
这个子程序是用来对任一方阵进行QR分解的，并且返回Q,R的值
）,q1（:
）,rq（:
real（kind=8）d,c0,h,he
allocate（u（n）,p（n）,w（n）,at（n,n）,c（n））
ma1（i,j）=ma（i,j）
q1（i,j）=0.0
else
q1（i,j）=1.0
endif
dor=1,n-1
he=he+ma1（i,r）**2
doi=r,n
if（ma1（r,r）>
h=c0**2-c0*ma1（r,r）
u（r）=ma1（r,r）-c0
u（r）=ma1（i,r）
at（i,j）=ma1（j,i）
callau（q1,u,c,n）
w（i）=c（i）
q1（i,j）=q1（i,j）-w（i）*u（j）/h
ma1（i,j）=ma1（i,j）-u（i）*p（j）
rq（i,j）=0.0
dor=1,n
rq（i,j）=rq（i,j）+ma1（i,r）*q1（r,j）
endsubroutineqr
subroutineopposite（matrix,u0,y,bo）!
这个子程序是用反幂法求出矩阵的模最小的特征值
real（kind=8）matrix（:
）,u0（:
）,y（:
）,bo
integern,i
real（kind=8）,allocatable:
）,us（:
）,x（:
real（kind=8）c,nn
bo=100
n=size（matrix,1）
allocate（u（n）,us（n,n）,x（n））
doi=1,n
u（i）=u0（i）
enddo
dowhile（1>
0）
c=bo
nn=sqrt（uu（u,u））
y（i）=u（i）/nn
calllu（matrix,us,y,x）!
调用对矩阵进行LU分解的子程序
u（i）=x（i）
enddo
bo=uu（y,u）
if（abs（bo-c）/abs（bo）<
=0.1e-3）exit
endsubroutineopposite
subroutinelu（matrix,us,b,x）!
这个子程序是用来对矩阵进行LU分解的，并且可以求解线性方程AX=B
）,b（:
l（:
real（kind=8）a,c,e,f
integeri,j,k,t
n=size（matrix,1）
allocate（l（n,n）,y（n））
a=0.0
us（1,j）=matrix（1,j）
doi=2,n
l（i,1）=matrix（i,1）/us（1,1）
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）
us（k,j）=matrix（k,j）-a
l（i,k）=（matrix（i,k）-c）/us（k,k）!
矩阵A的LU分解结束
dot=1,n-1
a=a+l（n,t）*us（t,n）
us（n,n）=matrix（n,n）-a
y
（1）=b
（1）
e=0.0
f=0.0
dot=1,i-1
e=e+l（i,t）*y（t）
y（i）=b（i）-e
e=0.0
x（n）=y（n）/us（n,n）
doi=n-1,1,-1
dot=i+1,n
f=f+us（i,t）*x（t）
x（i）=（y（i）-f）/us（i,i）!
求得解向量X（i）
endsubroutinelu
endmodule
programmain!
主程序开始的地方
uselinear2!
调用模块LINEAR2
implicitnone
integer,parameter:
nn=10
complex（kind=8）s1,s2,tr（nn）
real（kind=8）b（nn,nn）,b1（nn,nn）,d（2,2）,a（nn,nn）,mk（nn,nn）,a1（nn,nn）,q1（nn,nn）,ma1（nn,nn）,rq（nn,nn）
real（kind=8）u0（nn）,y（nn）,c（nn,nn）,bo,a11,a2
integeri,j,m
m=nn
doi=1,nn
u0（i）=i
doi=1,nn!
矩阵A的初始化
doj=1,nn
if（i/=j）then
b（i,j）=sin（0.5*i+0.2*j）
b（i,j）=1.5*cos（i+1.2*j）
calluptr（b,b1）!
对矩阵A进行拟上三角化，得矩阵A（N-1）
open（unit=30,file='
output/thicknessandf.txt'
）!
将程序运行的结果存储在一个TXT的文档里面
write（30,"
（'
矩阵A拟上三角化后所得的矩阵A（n-1）为：
'
）"
格式化输出A（N-1）
write（*,"
第'
i2,'
行前五个元素为'
5（1x,e20.12））"
）i,b1（i,1:
5）
行后五个元素为'
）i,b1（i,6:
10）
callqr（b1,ma1,q1,rq）!
对矩阵A（N-1）进行QR分解，并且输出Q,R,RQ等矩阵
矩阵A（n-1）进行QR分解后所得Q矩阵：
）i,q1（i,1:
）i,q1（i,6:
矩阵A（n-1）进行QR分解后所得R矩阵：
）i,ma1（i,1:
）i,ma1（i,6:
矩阵A（n-1）进行QR分解后所得的乘积矩阵RQ：
）i,rq（i,1:
）i,rq（i,6:
a（i,j）=b1（i,j）
dowhile（1.0>
0.0）
200if（abs（a（m,m-1））<
=0.1e-11）then!
对矩阵A（N-1）进行带双步位移的QR分解，求得其全部特征值
tr（m）=cmplx（a（m,m）,0.0）
矩阵的第'
个特征值'
'
tr（'
）=（'
e20.12,'
）'
）m,m,tr（m）
c（i,j）=sin（0.5*i+0.2*j）!
求原点平移后的矩阵A
c（i,j）=1.5*cos（i+1.2*j）-0.9*a（m,m）
callopposite（c,u0,y,bo）!
对实特征值运用带原点平移的反幂法，求得特征向量
A的相应于实特征值'
的特征向量为'
）a（m,m）
（e20.12）"
）y（i）
m=m-1
if（m==1）then
tr（m）=cmplx（a（1,1）,0.0）
c（i,j）=sin（0.5*i+0.2*j）!
callopposite（c,u0,y,bo）!
goto300
if（m==0）then
goto300
endif
if（m>
1）then
goto200
d（1,1）=a（m-1,m-1）
d（1,2）=a（m-1,m）
d（2,1）=a（m,m-1）
d（2,2）=a（m,m）
calltwo（d,s1,s2,a11,a2）!
求二阶子阵的两个特征值
if（m
升级会员 