计算方法Fortran版.docx
《计算方法Fortran版.docx》由会员分享，可在线阅读，更多相关《计算方法Fortran版.docx（49页珍藏版）》请在冰豆网上搜索。
计算方法Fortran版
一、 矩阵分解与线性方程组直接方法
1. LU分解解方程
moduleLU
!
----------------------------------------modulecoment
!
Version :
!
Codedby :
syz
!
Date :
!
-----------------------------------------------------
!
Description:
LU分解解方程
!
!
-----------------------------------------------------
contains
subroutinesolve（A,b,x,N）
implicitreal*8（a-z）
integer:
:
N
real*8:
:
A（N,N）,b（N）,x（N）
real*8:
:
L（N,N）,U（N,N）
real*8:
:
y（N）
calldoolittle（A,L,U,N）
call downtri（L,b,y,N）
calluptri（U,y,x,N）
endsubroutinesolve
subroutinedoolittle（A,L,U,N）
!
---------------------------------subroutine comment
!
Version :
!
Codedby :
syz
!
Date :
!
-----------------------------------------------------
!
Purpose :
LU分解之Doolittle函数
!
A=LU
!
-----------------------------------------------------
!
Input parameters :
!
1. A 方阵
!
2. N 阶数
!
Outputparameters :
!
1. L
!
2. U
!
Commonparameters :
!
!
----------------------------------------------------
!
PostScript:
!
1.
!
2.
!
----------------------------------------------------
implicitreal*8（a-z）
integer:
:
N,i,k,r
real*8:
:
A（N,N）,L（N,N）,U（N,N）
!
U的第一行
U（1,:
）=A（1,:
）
!
L的第一列
L（:
1）=a（:
1）/U（1,1）
dok=2,N
l（k,k）=1
doj=k,n
s=0
dom=1,k-1
s=s+l（k,m）*u（m,j）
enddo
u（k,j）=a（k,j）-s
enddo
doi=k+1,n
s=0
dom=1,k-1
s=s+l（i,m）*u（m,k）
enddo
l（i,k）=（a（i,k）-s）/u（k,k）
enddo
enddo
endsubroutinedoolittle
subroutineuptri（A,b,x,N）
!
---------------------------------subroutine comment
!
Version :
!
Codedby :
syz
!
Date :
2010-4-8
!
-----------------------------------------------------
!
Purpose :
上三角方程组的回带方法
!
Ax=b
!
-----------------------------------------------------
!
Input parameters :
!
1. A（N,N）系数矩阵
!
2. b（N）右向量
!
3. N方程维数
!
Outputparameters :
!
1. x 方程的根
!
2.
!
Commonparameters :
!
!
----------------------------------------------------
implicitreal*8（a-z）
integer:
:
i,j,k,N
real*8:
:
A（N,N）,b（N）,x（N）
x（N）=b（N）/A（N,N）
!
回带部分
doi=n-1,1,-1
x（i）=b（i）
doj=i+1,N
x（i）=x（i）-a（i,j）*x（j）
enddo
x（i）=x（i）/A（i,i）
enddo
endsubroutineuptri
subroutinedowntri（A,b,x,N）
!
---------------------------------subroutine comment
!
Version :
!
Codedby :
syz
!
Date :
2010-4-9
!
-----------------------------------------------------
!
Purpose :
下三角方程组的回带方法
!
Ax=b
!
-----------------------------------------------------
!
Input parameters :
!
1. A（N,N）系数矩阵
!
2. b（N）右向量
!
3. N方程维数
!
Outputparameters :
!
1. x 方程的根
!
2.
!
Commonparameters :
!
!
----------------------------------------------------
implicitreal*8（a-z）
integer:
:
i,j,N
real*8:
:
A（N,N）,b（N）,x（N）
x
（1）=b
（1）/a（1,1）
dok=2,N
x（k）=b（k）
doi=1,k-1
x（k）=x（k）-a（k,i）*x（i）
enddo
x（k）=x（k）/a（k,k）
enddo
endsubroutinedowntri
endmoduleLU
programmain
!
--------------------------------------programcomment
!
Version :
!
Codedby :
syz
!
Date :
2010-4-8
!
-----------------------------------------------------
!
Purpose :
LU分解计算线性方程组
!
Ax=b
!
-----------------------------------------------------
!
Inputdata files:
!
1. A,b
!
2.
!
Outputdatafiles :
!
1. x
!
2.
!
-----------------------------------------------------
useLU
integer,parameter:
:
N=4
real*8:
:
A（n,n）,L（N,N）,U（N,N）
real*8:
:
b（N）,x（N）
open（unit=11,file=''）
open（unit=12,file=''）
read（11,*）
read（11,*）（（A（i,j）,j=1,N）,i=1,N）
read（11,*）b
callsolve（A,b,x,N）
write（12,101）x
101format（T5,'LU分解计算线性方程组计算结果',LU分解之Crout
modulecrout
contains
subroutinesolve（A,L,U,N）
!
---------------------------------subroutine comment
!
Version :
!
Codedby :
syz
!
Date :
!
-----------------------------------------------------
!
Purpose :
LU之Crout分解
!
A=LU
!
-----------------------------------------------------
!
Input parameters :
!
1. A 方阵
!
2. N 阶数
!
Outputparameters :
!
1. L
!
2. U
!
Commonparameters :
!
!
----------------------------------------------------
!
PostScript:
!
1.
!
2.
!
----------------------------------------------------
implicitreal*8（a-z）
integer:
:
N,i,k,r
real*8:
:
A（N,N）,L（N,N）,U（N,N）
!
L第一列
L（:
1）=a（:
1）
!
U第一行
U（1,:
）=a（1,:
）/L（1,1）
dok=2,N
doi=k,n
s=0
dor=1,k-1
s=s+l（i,r）*u（r,k）
enddo
l（i,k）=a（i,k）-s
enddo
doj=k+1,n
s=0
dor=1,k-1
s=s+l（k,r）*u（r,j）
enddo
u（k,j）=（a（k,j）-s）/l（k,k）
enddo
u（k,k）=1
enddo
endsubroutinesolve
endmodulecrout
programmain
!
--------------------------------------programcomment
!
Version :
!
Codedby :
syz
!
Date :
2010-4-8
!
-----------------------------------------------------
!
Purpose :
Crot分解
!
!
-----------------------------------------------------
!
Inputdata files:
!
1. A,N
!
2.
!
Outputdatafiles :
!
1. L,U
!
2.
!
-----------------------------------------------------
usecrout
integer,parameter:
:
N=4
real*8:
:
A（n,n）,L（N,N）,U（N,N）
open（unit=11,file=''）
open（unit=12,file=''）
read（11,*）
read（11,*）（（A（i,j）,j=1,N）,i=1,N）
callsolve（A,L,U,N）
write（12,21）
21format（T10,'LU之Crout分解',/）
!
输出L矩阵
write（12,*）'L='
doi=1,N
write（12,22）L（i,:
）
enddo
22format
!
输出U矩阵
write（12,*）'U='
doi=1,N
write（12,22）U（i,:
）
enddo
23format
endprogrammain
3.LU分解之Doolittle
moduledoolittle
!
----------------------------------------modulecoment
!
Version :
!
Codedby :
syz
!
Date :
!
-----------------------------------------------------
!
Description:
LU分解之doolittle模块
!
!
-----------------------------------------------------
contains
subroutinesolve（A,L,U,N）
!
---------------------------------subroutine comment
!
Version :
!
Codedby :
syz
!
Date :
!
-----------------------------------------------------
!
Purpose :
LU分解之Doolittle函数
!
A=LU
!
-----------------------------------------------------
!
Input parameters :
!
1. A 方阵
!
2. N 阶数
!
Outputparameters :
!
1. L
!
2. U
!
Commonparameters :
!
!
----------------------------------------------------
!
PostScript:
!
1.
!
2.
!
----------------------------------------------------
implicitreal*8（a-z）
integer:
:
N,i,k,r
real*8:
:
A（N,N）,L（N,N）,U（N,N）
!
U的第一行
U（1,:
）=A（1,:
）
!
L的第一列
L（:
1）=a（:
1）/U（1,1）
dok=2,N
l（k,k）=1
doj=k,n
s=0
dom=1,k-1
s=s+l（k,m）*u（m,j）
enddo
u（k,j）=a（k,j）-s
enddo
doi=k+1,n
s=0
dom=1,k-1
s=s+l（i,m）*u（m,k）
enddo
l（i,k）=（a（i,k）-s）/u（k,k）
enddo
enddo
endsubroutinesolve
endmoduledoolittle
programmain
!
--------------------------------------programcomment
!
Version :
!
Codedby :
syz
!
Date :
2010-4-8
!
-----------------------------------------------------
!
Purpose :
Doolittle分解
!
!
-----------------------------------------------------
!
Inputdata files:
!
1. A,N
!
2.
!
Outputdatafiles :
!
1. L,U
!
2.
!
-----------------------------------------------------
usedoolittle
integer,parameter:
:
N=3
real*8:
:
A（n,n）,L（N,N）,U（N,N）
open（unit=11,file=''）
open（unit=12,file=''）
read（11,*）
read（11,*）（（A（i,j）,j=1,N）,i=1,N）
callsolve（A,L,U,N）
write（12,21）
21format（T10,'Doolittle分解',/）
!
输出L矩阵
write（12,*）'L='
doi=1,N
write（12,22）L（i,:
）
enddo
22format
!
输出U矩阵
write（12,*）'U='
doi=1,N
write（12,22）U（i,:
）
enddo
23format
endprogrammain
4.高斯消去法
modulegauss
!
----------------------------------------modulecoment
!
Version :
!
Codedby :
syz
!
Date :
2010-4-8
!
-----------------------------------------------------
!
Description:
高斯消去法模块
!
!
-----------------------------------------------------
!
Contains :
!
1. solve 方法函数
!
2.
!
-----------------------------------------------------
contains
subroutinesolve（A,b,x,N）
!
---------------------------------subroutine comment
!
Version :
!
Codedby :
syz
!
Date :
2010-4-8
!
-----------------------------------------------------
!
Purpose :
高斯消去法
!
Ax=b
!
-----------------------------------------------------
!
Input parameters :
!
1. A（N,N）系数矩阵
!
2. b（N）右向量
!
3. N方程维数
!
Outputparameters :
!
1. x 方程的根
!
2.
!
Commonparameters :
!
!
----------------------------------------------------
implicitreal*8（a-z）
integer:
:
i,k,N
real*8:
:
A（N,N）,b（N）,x（N）
real*8:
:
Aup（N,N）,bup（N）
!
Ab为增广矩阵 [Ab]
real*8:
:
Ab（N,N+1）
Ab（1:
N,1:
N）=A
Ab（:
N+1）=b
!
-------------------------------
!
这段是高斯消去法的核心部分
dok=1,N-1
doi=k+1,N
temp=Ab（i,k）/Ab（k,k）
Ab（i,:
）=Ab（i,:
）-temp*Ab（k,:
）
enddo
enddo
!
-----------------------------
!
经过上一步，Ab已经化为如下形式的矩阵
!
|* * * * #|
!
[Ab]=|0 * * * #|
!
|0 0 * * #|
!
|0 0 0 * #|
!
Aup（:
:
）=Ab（1:
N,1:
N）
bup（:
）=Ab（:
N+1）
!
调用用上三角方程组的回带方法
calluptri（Aup,bup,x,n）
endsubroutinesolve
subroutineuptri（A,b,x,N）
!
---------------------------------subroutine comment
!
Version :
!
Codedby :
syz
!
Date :
2010-4-8
!
-----------------------------------------------------
!
Purpose
升级会员 