计算方法Fortran版.docx

计算方法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