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