计算方法Fortran版Word下载.docx

计算方法Fortran版Word下载.docx

《计算方法Fortran版Word下载.docx》由会员分享，可在线阅读，更多相关《计算方法Fortran版Word下载.docx（49页珍藏版）》请在冰豆网上搜索。

integer:

:

N

real*8:

A（N,N）,b（N）,x（N）

L（N,N）,U（N,N）

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

Purpose 

LU分解之Doolittle函数

A=LU

Input 

parameters 

1. 

A 

方阵

2. 

N 

阶数

Outputparameters 

L

U

Commonparameters 

----------------------------------------------------

PostScript:

1.

2.

N,i,k,r

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

doi=k+1,n

s=s+l（i,m）*u（m,k）

l（i,k）=（a（i,k）-s）/u（k,k）

enddo

endsubroutinedoolittle

subroutineuptri（A,b,x,N）

2010-4-8

上三角方程组的回带方法

Ax=b

A（N,N）系数矩阵

b（N）右向量

3. 

N方程维数

x 

方程的根

i,j,k,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）

x（i）=x（i）/A（i,i）

endsubroutineuptri

subroutinedowntri（A,b,x,N）

2010-4-9

下三角方程组的回带方法

i,j,N

x

（1）=b

（1）/a（1,1）

x（k）=b（k）

doi=1,k-1

x（k）=x（k）-a（k,i）*x（i）

x（k）=x（k）/a（k,k）

endsubroutinedowntri

endmoduleLU

programmain

--------------------------------------programcomment

LU分解计算线性方程组

Inputdata 

files:

A,b

Outputdatafiles 

x

useLU

integer,parameter:

N=4

A（n,n）,L（N,N）,U（N,N）

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

subroutinesolve（A,L,U,N）

LU之Crout分解

L第一列

1）

U第一行

）=a（1,:

）/L（1,1）

doi=k,n

dor=1,k-1

s=s+l（i,r）*u（r,k）

l（i,k）=a（i,k）-s

doj=k+1,n

s=s+l（k,r）*u（r,j）

u（k,j）=（a（k,j）-s）/l（k,k）

u（k,k）=1

endmodulecrout

Crot分解

A,N

L,U

usecrout

callsolve（A,L,U,N）

write（12,21）

21format（T10,'

LU之Crout分解'

/）

输出L矩阵

write（12,*）'

L='

doi=1,N

write（12,22）L（i,:

22format

输出U矩阵

U='

write（12,22）U（i,:

23format

endprogrammain

3.LU分解之Doolittle

moduledoolittle

LU分解之doolittle模块

endmoduledoolittle

Doolittle分解

usedoolittle

N=3

Doolittle分解'

4.高斯消去法

modulegauss

高斯消去法模块

Contains 

solve 

方法函数

contains

高斯消去法

i,k,N

Aup（N,N）,bup（N）

Ab为增广矩阵 

[Ab]

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,:

-----------------------------

经过上一步，Ab已经化为如下形式的矩阵

|* 

* 

#|

[Ab]=|0 

|0 

0 

Aup（:

:

）=Ab（1:

N）

bup（:

）=Ab（:

N+1）

调用用上三角方程组的回带方法

calluptri（Aup,bup,x,n）