计算方法Fortran版.docx
《计算方法Fortran版.docx》由会员分享，可在线阅读，更多相关《计算方法Fortran版.docx（19页珍藏版）》请在冰豆网上搜索。
计算方法Fortran版
-CAL-FENGHAI.-（YICAI）-CompanyOne1-CAL-本页仅作为文档封面，使用请直接删除
计算方法Fortran版（总20页）
第一章：
线性方程组的解法
一、矩阵分解与线性方程组直接方法
1.LU分解解方程
moduleLU
!
----------------------------------------modulecoment!
Version:
!
Codedby:
syz!
Date:
!
-----------------------------------------------------!
Description:
LU分解解方程!
!
-----------------------------------------------------
contains?
subroutinesolve（A,b,x,N）
implicitreal*8（a-z）integer:
:
Nreal*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）calldowntri（L,b,y,N）calluptri（U,y,x,N）
endsubroutinesolve
subroutinedoolittle（A,L,U,N）!
---------------------------------subroutinecomment!
Version:
!
Codedby:
syz!
Date:
!
-----------------------------------------------------!
Purpose:
LU分解之Doolittle函数!
A=LU!
-----------------------------------------------------!
Inputparameters:
!
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,Nl（k,k）=1doj=k,ns=0dom=1,k-1s=s+l（k,m）*u（m,j）enddou（k,j）=a（k,j）-senddodoi=k+1,ns=0dom=1,k-1s=s+l（i,m）*u（m,k）enddol（i,k）=（a（i,k）-s）/u（k,k）enddoenddo
endsubroutinedoolittle
subroutineuptri（A,b,x,N）!
---------------------------------subroutinecomment!
Version:
!
Codedby:
syz!
Date:
2010-4-8!
-----------------------------------------------------!
Purpose:
上三角方程组的回带方法!
Ax=b!
-----------------------------------------------------!
Inputparameters:
!
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,-1x（i）=b（i）doj=i+1,Nx（i）=x（i）-a（i,j）*x（j）enddox（i）=x（i）/A（i,i）
enddo
endsubroutineuptri
subroutinedowntri（A,b,x,N）!
---------------------------------subroutinecomment!
Version:
!
Codedby:
syz!
Date:
2010-4-9!
-----------------------------------------------------!
Purpose:
下三角方程组的回带方法!
Ax=b!
-----------------------------------------------------!
Inputparameters:
!
1.A（N,N）系数矩阵!
2.b（N）右向量!
3.N方程维数!
Outputparameters:
!
1.x方程的根!
2.!
Commonparameters:
!
!
----------------------------------------------------
implicitreal*8（a-z）integer:
:
i,j,Nreal*8:
:
A（N,N）,b（N）,x（N）
x
（1）=b
（1）/a（1,1）
dok=2,Nx（k）=b（k）doi=1,k-1x（k）=x（k）-a（k,i）*x（i）enddox（k）=x（k）/a（k,k）
enddo
endsubroutinedowntri
endmoduleLU
?
programmain!
--------------------------------------programcomment!
Version:
!
Codedby:
syz!
Date:
2010-4-8!
-----------------------------------------------------!
Purpose:
LU分解计算线性方程组!
Ax=b!
-----------------------------------------------------!
Inputdatafiles:
!
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）!
---------------------------------subroutinecomment!
Version:
!
Codedby:
syz!
Date:
!
-----------------------------------------------------!
Purpose:
LU之Crout分解!
A=LU!
-----------------------------------------------------!
Inputparameters:
!
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,ns=0dor=1,k-1s=s+l（i,r）*u（r,k）enddol（i,k）=a（i,k）-senddodoj=k+1,ns=0dor=1,k-1s=s+l（k,r）*u（r,j）enddou（k,j）=（a（k,j）-s）/l（k,k）enddou（k,k）=1enddo
endsubroutinesolve
endmodulecrout
?
programmain!
--------------------------------------programcomment!
Version:
!
Codedby:
syz!
Date:
2010-4-8!
-----------------------------------------------------!
Purpose:
Crot分解!
!
-----------------------------------------------------!
Inputdatafiles:
!
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,:
）enddo22format
!
输出U矩阵write（12,*）'U='doi=1,N?
write（12,22）U（i,:
）enddo23format
endprogrammain
3.LU分解之Doolittle
moduledoolittle
!
----------------------------------------modulecoment!
Version:
!
Codedby:
syz!
Date:
!
-----------------------------------------------------!
Description:
LU分解之doolittle模块!
!
-----------------------------------------------------
contains?
subroutinesolve（A,L,U,N）!
---------------------------------subroutinecomment!
Version:
!
Codedby:
syz!
Date:
!
-----------------------------------------------------!
Purpose:
LU分解之Doolittle函数!
A=LU!
-----------------------------------------------------!
Inputparameters:
!
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,Nl（k,k）=1doj=k,ns=0dom=1,k-1s=s+l（k,m）*u（m,j）enddou（k,j）=a（k,j）-senddodoi=k+1,ns=0dom=1,k-1s=s+l（i,m）*u（m,k）enddol（i,k）=（a（i,k）-s）/u（k,k）enddoenddo
endsubroutinesolve
endmoduledoolittle
?
programmain!
--------------------------------------programcomment!
Version:
!
Codedby:
syz!
Date:
2010-4-8!
-----------------------------------------------------!
Purpose:
Doolittle分解!
!
-----------------------------------------------------!
Inputdatafiles:
!
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,:
）enddo22format
!
输出U矩阵write（12,*）'U='doi=1,N?
write（12,22）U（i,:
）enddo23format
endprogrammain
4.高斯消去法
modulegauss!
----------------------------------------modulecoment!
Version:
!
Codedby:
syz!
Date:
2010-4-8!
-----------------------------------------------------!
Description:
高斯消去法模块!
!
-----------------------------------------------------!
Contains:
!
1.solve方法函数!
2.!
-----------------------------------------------------
contains
subroutinesolve（A,b,x,N）!
---------------------------------subroutinecomment!
Version:
!
Codedby:
syz!
Date:
2010-4-8!
-----------------------------------------------------!
Purpose:
高斯消去法!
Ax=b!
-----------------------------------------------------!
Inputparameters:
!
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,Ntemp=Ab（i,k）/Ab（k,k）Ab（i,:
）=Ab（i,:
）-temp*Ab（k,:
）enddo
enddo
!
-----------------------------!
经过上一步，Ab已经化为如下形式的矩阵!
|****#|!
[Ab]=|0***#|!
|00**#|!
|000*#|!
Aup（:
:
）=Ab（1:
N,1:
N）
bup（:
）=Ab（:
N+1）
!
调用用上三角方程组的回带方法calluptri（Aup,bup,x,n）
endsubroutinesolve
?
subroutineuptri（A,b,x,N）!
---------------------------------subroutinecomment!
Version:
!
Codedby:
syz!
Date:
2010-4-8!
-----------------------------------------------------!
Purpose:
上三角方程组的回带方法!
Ax=b!
-----------------------------------------------------!
Inputparameters:
!
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（N）=b（N）/A（N,N）
!
回带部分doi=n-1,1,-1x（i）=b（i）doj=i+1,Nx（i）=x（i）-a（i,j）*x（j）enddox（i）=x（i）/A（i,i）
enddo
endsubroutineuptri
endmodulegauss
programmain!
--------------------------------------programcomment!
Version:
!
Codedby:
syz!
Date:
2010-4-8!
-----------------------------------------------------!
Purpose:
高斯消去法!
!
-----------------------------------------------------!
Inputdatafiles:
!
1.输入方程系数!
2.!
Outputdatafiles:
!
1.计算结果!
2.!
-----------------------------------------------------!
PostScript:
!
1.需要准备输入数据!
!
-----------------------------------------------------usegauss
implicitreal*8（a-z）
integer,parameter:
:
N=4
integer:
:
i,jreal*8:
:
A（N,N）,b（N）,x（N）
?
open（unit=11,file=''）open（unit=12,file=''）
read（11,*）!
读入A矩阵
read（11,*）（（A（i,j）,j=1,N）,i=1,N）!
读入B向量read（11,*）b
callsolve（A,b,x,N）
write（12,101）x101format（T5,'高斯消去法计算结果',/,T4,'x=',4（/）
endprogrammain
5.列主元消去法
modulem_gauss!
----------------------------------------modulecoment!
Version:
!
Codedby:
syz!
Date:
2010-4-8!
-----------------------------------------------------!
Description:
高斯列主元消去法模块!
!
-----------------------------------------------------!
Contains:
!
1.solve方法函数!
2.!
-----------------------------------------------------
contains
subroutinesolve（A,b,x,N）!
---------------------------------subroutinecomment!
Version:
!
Codedby:
syz!
Date:
2010-4-8!
-----------------------------------------------------!
Purpose:
高斯列主元消去法!
Ax=b!
-----------------------------------------------------!
Inputp
升级会员 