北航数值分析大作业2.docx
《北航数值分析大作业2.docx》由会员分享,可在线阅读,更多相关《北航数值分析大作业2.docx(30页珍藏版)》请在冰豆网上搜索。
北航数值分析大作业2
《数值分析B》
第二次数值分析大作业
院系:
04能源与动力工程学院
姓名:
王开逍
学号:
SY1104207
一.算法设计方案:
1、对矩阵A进行拟上三角化,得矩阵A(N-1).
2、对矩阵A(N-1)运用带双步位移的QR分解方法,求得其全部特征值,该特征值就是矩阵A的全部特征值.
3、对矩阵A(N-1)进行QR分解,求得矩阵Q,R,RQ.
4、对于矩阵的实特征值,通过带原点平移的反幂法,求得其特征向量.
二.源程序:
该计算程序使用FORTRAN语言编写
modulelinear2
implicitnone
contains
subroutineau(matrix,u,c,n)!
这个子程序是用来计算矩阵和向量相乘的,最后返回一个向量值
real(kind=8)matrix(:
:
),u(:
),c(:
)
integern,i,k
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
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
subroutineuptr(ma,ma1)!
这个子程序是用来对矩阵A进行拟上三角化的
real(kind=8)ma(:
:
),ma1(:
:
)
integern
real(kind=8),allocatable:
:
u(:
),p(:
),q(:
),w(:
),wu(:
:
),up(:
:
),at(:
:
),c(:
)
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))
doi=1,n
doj=1,n
ma1(i,j)=0.0
enddo
enddo!
ma1的零初始化
dor=1,n-2
he=0.0
doi=r+1,n
he=he+ma(i,r)**2
enddo
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)
enddo
ma1(r+1,r)=c0
he=0.0
doi=r+2,n
he=he+ma(i,r)**2
enddo
if(he==0)then
cycle
else
h=c0**2-c0*ma(r+1,r)
doi=1,r
u(i)=0.0
enddo
u(r+1)=ma(r+1,r)-c0
doi=r+2,n
u(i)=ma(i,r)
enddo
doi=1,n
doj=1,n
at(i,j)=ma(j,i)
enddo
enddo
callau(at,u,c,n)
doi=1,n
p(i)=c(i)/h
enddo
callau(ma,u,c,n)
doi=1,n
q(i)=c(i)/h
enddo
t=uu(p,u)/h
doi=1,n
w(i)=q(i)-t*u(i)
enddo
doi=1,n
doj=1,n
ma(i,j)=ma(i,j)-w(i)*u(j)-u(i)*p(j)
enddo
enddo
endif
enddo
doi=1,n
ma1(i,n-1)=ma(i,n-1)
enddo
doi=1,n
ma1(i,n)=ma(i,n)
enddo
endsubroutineuptr
subroutinetwo(d,s1,s2,a1,a2)!
这个子程序是用来求二阶子阵的两个特征值的
real(kind=8)d(:
:
)
complex(kind=8)s1,s2
real(kind=8)a1,a2,b1,b2
integern
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
else
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
endif
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)
enddo
if(i/=j)then
mk(i,j)=c-s*ak(i,j)
else
mk(i,j)=c-s*ak(i,j)+t0
endif
enddo
enddo
callmkak(ak,mk,ak1,m)
endsubroutineak2
subroutinemkak(ak,mk,ak1,m)!
!
用来对中间矩阵M(K)进行QR分解的,并且求出A(K+1)
real(kind=8)mk(:
:
),ak(:
:
),ak1(:
:
)
integerm
integeri,j,r
real(kind=8)d,c0,h,t,he
real(kind=8),allocatable:
:
b(:
:
),c1(:
:
),u(:
),v(:
),bt(:
:
),c1t(:
:
),p(:
),q(:
),w(:
),c(:
)
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))
doi=1,m
doj=1,m
b(i,j)=mk(i,j)
c1(i,j)=ak(i,j)
enddo
enddo
dor=1,m-1
he=0.0
doi=r+1,m
he=he+b(i,r)**2
enddo
if(he==0)then
cycle
else
he=0.0
doi=r,m
he=he+b(i,r)**2
enddo
d=sqrt(he)
if(b(r,r)>0.0)then
c0=0.0-d
else
c0=d
endif
h=c0**2-c0*b(r,r)
doi=1,m
u(i)=0.0
enddo
u(r)=b(r,r)-c0
doi=r+1,m
u(i)=b(i,r)
enddo
doi=1,m
doj=1,m
bt(i,j)=b(j,i)
c1t(i,j)=c1(j,i)
enddo
enddo
callau(bt,u,c,m)
doi=1,m
v(i)=c(i)/h
enddo
doi=1,m
doj=1,m
b(i,j)=b(i,j)-u(i)*v(j)
enddo
enddo
callau(c1t,u,c,m)
doi=1,m
p(i)=c(i)/h
enddo
callau(c1,u,c,m)
doi=1,m
q(i)=c(i)/h
enddo
t=uu(p,u)/h
doi=1,m
w(i)=q(i)-t*u(i)
enddo
doi=1,m
doj=1,m
c1(i,j)=c1(i,j)-w(i)*u(j)-u(i)*p(j)
enddo
enddo
endif
enddo
doi=1,m
doj=1,m
ak1(i,j)=c1(i,j)
enddo
enddo
endsubroutinemkak
subroutineqr(ma,ma1,q1,rq)!
这个子程序是用来对任一方阵进行QR分解的,并且返回Q,R的值
real(kind=8)ma(:
:
),ma1(:
:
),q1(:
:
),rq(:
:
)
integern
real(kind=8),allocatable:
:
u(:
),p(:
),w(:
),at(:
:
),c(:
)
real(kind=8)d,c0,h,he
integeri,j,r
n=size(ma,1)
allocate(u(n),p(n),w(n),at(n,n),c(n))
doi=1,n
doj=1,n
ma1(i,j)=ma(i,j)
enddo
enddo
doi=1,n
doj=1,n
if(i/=j)then
q1(i,j)=0.0
else
q1(i,j)=1.0
endif
enddo
enddo
dor=1,n-1
he=0.0
doi=r+1,n
he=he+ma1(i,r)**2
enddo
if(he==0)then
cycle
else
he=0.0
doi=r,n
he=he+ma1(i,r)**2
enddo
d=sqrt(he)
if(ma1(r,r)>0.0)then
c0=0.0-d
else
c0=d
endif
h=c0**2-c0*ma1(r,r)
doi=1,n
u(i)=0.0
enddo
u(r)=ma1(r,r)-c0
doi=r+1,n
u(r)=ma1(i,r)
enddo
doi=1,n
doj=1,n
at(i,j)=ma1(j,i)
enddo
enddo
callau(q1,u,c,n)
doi=1,n
w(i)=c(i)
enddo
doi=1,n
doj=1,n
q1(i,j)=q1(i,j)-w(i)*u(j)/h
enddo
enddo
callau(at,u,c,n)
doi=1,n
p(i)=c(i)/h
enddo
doi=1,n
doj=1,n
ma1(i,j)=ma1(i,j)-u(i)*p(j)
enddo
enddo
endif
enddo
doi=1,n
doj=1,n
rq(i,j)=0.0
dor=1,n
rq(i,j)=rq(i,j)+ma1(i,r)*q1(r,j)
enddo
enddo
enddo
endsubroutineqr
subroutineopposite(matrix,u0,y,bo)!
这个子程序是用反幂法求出矩阵的模最小的特征值
real(kind=8)matrix(:
:
),u0(:
),y(:
),bo
integern,i
real(kind=8),allocatable:
:
u(:
),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))
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)
if(abs(bo-c)/abs(bo)<=0.1e-3)exit
enddo
endsubroutineopposite
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
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
enddo
doi=1,nn!
!
矩阵A的初始化
doj=1,nn
if(i/=j)then
b(i,j)=sin(0.5*i+0.2*j)
else
b(i,j)=1.5*cos(i+1.2*j)
endif
enddo
enddo
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(*,"('矩阵A拟上三角化后所得的矩阵A(n-1)为:
')")
doi=1,m
write(30,"('第',i2,'行前五个元素为',5(1x,e20.12))")i,b1(i,1:
5)
write(30,"('第',i2,'行后五个元素为',5(1x,e20.12))")i,b1(i,6:
10)
write(*,"('第',i2,'行前五个元素为',5(1x,e20.12))")i,b1(i,1:
5)
write(*,"('第',i2,'行后五个元素为',5(1x,e20.12))")i,b1(i,6:
10)
enddo
callqr(b1,ma1,q1,rq)!
对矩阵A(N-1)进行QR分解,并且输出Q,R,RQ等矩阵
write(30,"('矩阵A(n-1)进行QR分解后所得Q矩阵:
')")
write(*,"('矩阵A(n-1)进行QR分解后所得Q矩阵:
')")
doi=1,m
write(30,"('第',i2,'行前五个元素为',5(1x,e20.12))")i,q1(i,1:
5)
write(30,"('第',i2,'行后五个元素为',5(1x,e20.12))")i,q1(i,6:
10)
write(*,"('第',i2,'行前五个元素为',5(1x,e20.12))")i,q1(i,1:
5)
write(*,"('第',i2,'行后五个元素为',5(1x,e20.12))")i,q1(i,6:
10)
enddo
write(30,"('矩阵A(n-1)进行QR分解后所得R矩阵:
')")
write(*,"('矩阵A(n-1)进行QR分解后所得R矩阵:
')")
doi=1,m
write(30,"('第',i2,'行前五个元素为',5(1x,e20.12))")i,ma1(i,1:
5)
write(30,"('第',i2,'行后五个元素为',5(1x,e20.12))")i,ma1(i,6:
10)
write(*,"('第',i2,'行前五个元素为',5(1x,e20.12))")i,ma1(i,1:
5)
write(*,"('第',i2,'行后五个元素为',5(1x,e20.12))")i,ma1(i,6:
10)
enddo
write(30,"('矩阵A(n-1)进行QR分解后所得的乘积矩阵RQ:
')")
write(*,"('矩阵A(n-1)进行QR分解后所得的乘积矩阵RQ:
')")
doi=1,m
write(30,"('第',i2,'行前五个元素为',5(1x,e20.12))")i,rq(i,1:
5)
write(30,"('第',i2,'行后五个元素为',5(1x,e20.12))")i,rq(i,6:
10)
write(*,"('第',i2,'行前五个元素为',5(1x,e20.12))")i,rq(i,1:
5)
write(*,"('第',i2,'行后五个元素为',5(1x,e20.12))")i,rq(i,6:
10)
enddo
doi=1,m
doj=1,m
a(i,j)=b1(i,j)
enddo
enddo
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)
write(30,"('矩阵的第',i2,'个特征值','tr(',i2,')=(',e20.12,',',e20.12,')')")m,m,tr(m)
write(*,"('矩阵的第',i2,'个特征值','tr(',i2,')=(',e20.12,',',e20.12,')')")m,m,tr(m)
doi=1,nn
doj=1,nn
if(i/=j)then
c(i,j)=sin(0.5*i+0.2*j)!
求原点平移后的矩阵A
else
c(i,j)=1.5*cos(i+1.2*j)-0.9*a(m,m)
endif
enddo
enddo
callopposite(c,u0,y,bo)!
对实特征值运用带原点平移的反幂法,求得特征向量
write(30,"('A的相应于实特征值',e20.12,'的特征向量为')")a(m,m)
write(*,"('A的相应于实特征值',e20.12,'的特征向量为')")a(m,m)
doi=1,nn
write(30,"(e20.12)")y(i)
write(*,"(e20.12)")y(i)
enddo
m=m-1
if(m==1)then
tr(m)=cmplx(a(1,1),0.0)
write(30,"('矩阵的第',i2,'个特征值','tr(',i2,')=(',e20.12,',',e20.12,')')")m,m,tr(m)
write(*,"('矩阵的第',i2,'个特征值','tr(',i2,')=(',e20.12,',',e20.12,')')")m,m,tr(m)
doi=1,nn
doj=1,nn
if(i/=j)then
c(i,j)=sin(0.5*i+0.2*j)!
求原点平移后的矩阵A
else
c(i,j)=1.5*cos(i+1.2*j)-0.9*a(m,m)
endif
enddo
enddo
callopposite(c,u0,y,bo)!
对实特征值运用带原点平移的反幂法,求得特征向量
write(30,"('A的相应于实特征值',e20.12,'的特征向量为')")a(m,m)
write(*,"('A的相应于实特征值',e20.12,'的特征向量为')")a(m,m)
doi=1,nn
write(30,"(e20.12)")y(i)
write(*,"(e20.12)")y(i)
enddo
goto300
endif
if(m==0)then
goto300
endif
if(m>1)then
goto200
endif
else
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