北航数值分析大作业2Word文档格式.docx
- 文档编号:21383700
- 上传时间:2023-01-30
- 格式:DOCX
- 页数:30
- 大小:93.06KB
北航数值分析大作业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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 作业