《计算方法》Fortran版

1、第一章：线性方程组的解法一 、 矩阵分解与线性方程组直接方法 1.  LU分解解方程module LU!-module coment!  Version     :  V1.0    !  Coded by    :  syz !  Date        :  !-!  Description :

2、  LU 分解解方程!    !-contains subroutine solve(A,b,x,N)implicit real*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) call doolittle(A,L,U,N)   call  downtri(L,b,y,N)  call uptri(U,y,x,N)end subroutine solvesubroutine

3、doolittle(A,L,U,N)!-subroutine  comment!  Version   :  V1.0    !  Coded by  :  syz !  Date      :  !-!  Purpose   :  LU分解之Doolittle函数!       

4、;       A=LU!-!  Input  parameters  :!       1.    A  方阵!       2.    N  阶数!  Output parameters  :!       1. &

5、#160; L!       2.   U!  Common parameters  :!-!  Post Script :!       1.!       2.!-implicit real*8(a-z)integer:N,i,k,rreal*8:A(N,N),L(N,N),U(N,N)!U的第一行U(1,:)=A(1,:)!L的第一列L(:,1)=a(:,1)/

6、U(1,1)do k=2,N       l(k,k)=1      do j=k,n       s=0       do m=1,k-1        s=s+l(k,m)*u(m,j)       end do  

7、     u(k,j)=a(k,j)-s   end do         do i=k+1,n     s=0     do m=1,k-1      s=s+l(i,m)*u(m,k)     end do     l(i,k)=(a(i,k)-

8、s)/u(k,k)          end do    end doend subroutine doolittlesubroutine uptri(A,b,x,N)!-subroutine  comment!  Version   :  V1.0    !  Coded by  :  syz !  Date 

9、60;    :  2010-4-8!-!  Purpose   :  上三角方程组的回带方法!                 Ax=b!-!  Input  parameters  :!       1.   A(N,N)系数矩阵! &

10、#160;     2.   b(N)右向量!       3.   N方程维数!  Output parameters  :!       1.  x  方程的根!       2.!  Common parameters  :!-implicit real*8(a-z)int

11、eger:i,j,k,Nreal*8:A(N,N),b(N),x(N)x(N)=b(N)/A(N,N)!回带部分do i=n-1,1,-1       x(i)=b(i)   do j=i+1,N    x(i)=x(i)-a(i,j)*x(j)   end do    x(i)=x(i)/A(i,i)end doend subroutine uptrisubroutine downtri(A,b,x,N)!-subroutine

12、60; comment!  Version   :  V1.0    !  Coded by  :  syz !  Date      :  2010-4-9!-!  Purpose   :  下三角方程组的回带方法!            

13、0;    Ax=b!-!  Input  parameters  :!       1.   A(N,N)系数矩阵!       2.   b(N)右向量!       3.   N方程维数!  Output parameters  :!   &#

14、160;   1.  x  方程的根!       2.!  Common parameters  :!-implicit real*8(a-z)integer:i,j,Nreal*8:A(N,N),b(N),x(N)x(1)=b(1)/a(1,1)do k=2,N   x(k)=b(k)   do i=1,k-1      x(k)=x(k)-a(k,i)*x(i) 

15、0; end do   x(k)=x(k)/a(k,k)end doend subroutine downtriend module LU program main!-program comment!  Version   :  V1.0    !  Coded by  :  syz !  Date      :  2010-4-8!-!  Purpose 

16、  :   LU分解计算线性方程组!              Ax=b!-!  In put data  files :!       1.    A,b!       2.!  Output data files  :!  &

17、#160;    1.    x!       2.!-use LUinteger,parameter:N=4real*8:A(n,n),L(N,N),U(N,N)real*8:b(N),x(N)open(unit=11,file='fin.txt')open(unit=12,file='fout.txt')read(11,*)read(11,*)(A(i,j),j=1,N),i=1,N)read(11,*)bcall solve(A,b,x,

18、N)write(12,101)x101 format(T5,'LU分解计算线性方程组计算结果',/,4(/,F10.6)end program main2. LU分解之Croutmodule croutcontains subroutine solve(A,L,U,N)!-subroutine  comment!  Version   :  V1.0    !  Coded by  :  syz !  Date 

19、0;    :  !-!  Purpose   :  LU之Crout分解!              A=LU!-!  Input  parameters  :!       1.    A  方阵!    

20、   2.    N  阶数!  Output parameters  :!       1.   L!       2.   U!  Common parameters  :!-!  Post Script :!       1.!   

21、    2.!-implicit real*8(a-z)integer:N,i,k,rreal*8:A(N,N),L(N,N),U(N,N)!L 第一列L(:,1)=a(:,1)!U 第一行U(1,:)=a(1,:)/L(1,1) do k=2,N   do i=k,n       s=0       do r=1,k-1        s=s+l

22、(i,r)*u(r,k)       end do       l(i,k)=a(i,k)-s   end do         do j=k+1,n     s=0     do r=1,k-1      s=s+l(k,r)*u(r,j)

23、0;    end do     u(k,j)=(a(k,j)-s)/l(k,k)          end do   u(k,k)=1   end doend subroutine solveend module crout program main!-program comment!  Version   :  V1.0 &

24、#160;  !  Coded by  :  syz !  Date      :  2010-4-8!-!  Purpose   :   Crot 分解!    !-!  In put data  files :!       1.    A,N! 

25、0;     2.!  Output data files  :!       1.    L,U!       2.!-use croutinteger,parameter:N=4real*8:A(n,n),L(N,N),U(N,N)open(unit=11,file='fin.txt')open(unit=12,file='fout.txt')read

26、(11,*)read(11,*)(A(i,j),j=1,N),i=1,N)call solve(A,L,U,N)write(12,21)21 format(T10,'LU之Crout分解',/)!输出L矩阵write(12,*)'L='do i=1,N write(12,22)L(i,:)end do22 format(4F10.6)!输出U矩阵write(12,*)'U='do i=1,N write(12,22)U(i,:)end do23 format(4F10.6)end program main3. LU分解之Doo

27、littlemodule doolittle!-module coment!  Version     :  V1.0    !  Coded by    :  syz !  Date        :  !-!  Description :  LU 分解之doolittle模块!   &

28、#160;!-contains subroutine solve(A,L,U,N)!-subroutine  comment!  Version   :  V1.0    !  Coded by  :  syz !  Date      :  !-!  Purpose   :  LU分解之Doolittle函数!  &

29、#160;           A=LU!-!  Input  parameters  :!       1.    A  方阵!       2.    N  阶数!  Output parameters  :!   

30、    1.   L!       2.   U!  Common parameters  :!-!  Post Script :!       1.!       2.!-implicit real*8(a-z)integer:N,i,k,rreal*8:A(N,N),L(N,N),U(N,N)!U的第一行U(1,:

31、)=A(1,:)!L的第一列L(:,1)=a(:,1)/U(1,1) do k=2,N       l(k,k)=1      do j=k,n       s=0       do m=1,k-1        s=s+l(k,m)*u(m,j)   

32、60;   end do       u(k,j)=a(k,j)-s   end do         do i=k+1,n     s=0     do m=1,k-1      s=s+l(i,m)*u(m,k)     end do 

33、    l(i,k)=(a(i,k)-s)/u(k,k)          end do    end doend subroutine solveend module doolittle program main!-program comment!  Version   :  V1.0    !  Coded by  :

34、60; syz !  Date      :  2010-4-8!-!  Purpose   :   Doolittle 分解!    !-!  In put data  files :!       1.    A,N!       2.!  Ou

35、tput data files  :!       1.    L,U!       2.!-use doolittleinteger,parameter:N=3real*8:A(n,n),L(N,N),U(N,N)open(unit=11,file='fin.txt')open(unit=12,file='fout.txt')read(11,*)read(11,*)(A(i,j),j=1,N),i=1,

36、N)call solve(A,L,U,N)write(12,21)21 format(T10,'Doolittle分解',/)!输出L矩阵write(12,*)'L='do i=1,N write(12,22)L(i,:)end do22 format(3F10.6)!输出U矩阵write(12,*)'U='do i=1,N write(12,22)U(i,:)end do23 format(3F10.6)end program main4.高斯消去法module gauss!-module coment!  Ver

37、sion     :  V1.0    !  Coded by    :  syz !  Date        :  2010-4-8!-!  Description : 高斯消去法模块!    !-!  Contains    :!   

38、;   1.   solve  方法函数!      2.!-containssubroutine solve(A,b,x,N)!-subroutine  comment!  Version   :  V1.0    !  Coded by  :  syz !  Date      :  2010-

39、4-8!-!  Purpose   :  高斯消去法!                 Ax=b!-!  Input  parameters  :!       1.   A(N,N)系数矩阵!       2. 

40、60; b(N)右向量!       3.   N方程维数!  Output parameters  :!       1.  x  方程的根!       2.!  Common parameters  :!-implicit real*8(a-z)integer:i,k,Nreal*8:A(N,N),b(N),x(N)real*8:A

41、up(N,N),bup(N)!Ab为增广矩阵  Abreal*8:Ab(N,N+1)Ab(1:N,1:N)=AAb(:,N+1)=b!-!  这段是 高斯消去法的核心部分do k=1,N-1   do i=k+1,N        temp=Ab(i,k)/Ab(k,k)          Ab(i,:)=Ab(i,:)-temp*Ab(k,:)    &#

42、160; end doend do!-! 经过上一步，Ab已经化为如下形式的矩阵!            | *  *  *  *  # |!     A b= | 0  *  *  *  # |!            | 0  0  *

43、  *  # |!            | 0  0  0  *  # |!Aup(:,:)=Ab(1:N,1:N)bup(:)=Ab(:,N+1)!调用用上三角方程组的回带方法call uptri(Aup,bup,x,n)end subroutine solve subroutine uptri(A,b,x,N)!-subroutine  comment!  Version 

44、0; :  V1.0    !  Coded by  :  syz !  Date      :  2010-4-8!-!  Purpose   :  上三角方程组的回带方法!                 Ax=b!-!  I

45、nput  parameters  :!       1.   A(N,N)系数矩阵!       2.   b(N)右向量!       3.   N方程维数!  Output parameters  :!       1.  x  方程

46、的根!       2.!  Common parameters  :!-implicit real*8(a-z)integer:i,j,Nreal*8:A(N,N),b(N),x(N)x(N)=b(N)/A(N,N)!回带部分do i=n-1,1,-1       x(i)=b(i)   do j=i+1,N    x(i)=x(i)-a(i,j)*x(j)   end do 

47、;   x(i)=x(i)/A(i,i)end doend subroutine uptriend module gaussprogram main!-program comment!  Version   :  V1.0    !  Coded by  :  syz !  Date      :  2010-4-8!-!  Purpose   : 

48、; 高斯消去法!    !-!  In put data  files :!       1.  fin.txt  输入方程系数!       2.!  Output data files  :!       1. fout.txt  计算结果!     &#

49、160; 2.!-!  Post Script :!       1.    需要准备输入数据!     !-use gaussimplicit real*8(a-z)integer,parameter: N=4integer:i,jreal*8:A(N,N),b(N),x(N) open(unit=11,file='fin.txt')open(unit=12,file='fout.txt')read(1

50、1,*)!读入A矩阵read(11,*)(A(i,j),j=1,N),i=1,N)!读入B向量read(11,*) bcall solve(A,b,x,N)write(12,101)x101 format(T5,'高斯消去法计算结果',/,T4,'x=',4(/F12.8)end program main5.列主元消去法module m_gauss!-module coment!  Version     :  V1.0    !  Coded by&#

51、160;   :  syz !  Date        :  2010-4-8!-!  Description : 高斯列主元消去法模块!    !-!  Contains    :!      1.   solve  方法函数!      2.!

52、-containssubroutine solve(A,b,x,N)!-subroutine  comment!  Version   :  V1.0    !  Coded by  :  syz !  Date      :  2010-4-8!-!  Purpose   :  高斯列主元消去法!     

53、            Ax=b!-!  Input  parameters  :!       1.   A(N,N)系数矩阵!       2.   b(N)右向量!       3.   N方程维数! 

54、Output parameters  :!       1.  x  方程的根!       2.!  Common parameters  :!-implicit real*8(a-z)integer:i,k,Ninteger:id_max  !主元素标号real*8:A(N,N),b(N),x(N)real*8:Aup(N,N),bup(N)!Ab为增广矩阵  Abreal*8:Ab(N,N+1)rea

55、l*8:vtemp1(N+1),vtemp2(N+1)Ab(1:N,1:N)=AAb(:,N+1)=b!#!  这段是 列主元消去法的核心部分do k=1,N-1    elmax=dabs(Ab(k,k)    id_max=k        !这段为查找主元素     !这段程序的主要目的不是为了赋值最大元素给elmax，而是为了找出最大元素对应的标号  do i=k+1,n 

56、60;    if (dabs(Ab(i,k)>elmax) then         elmax=Ab(i,k)         id_max=i      end if              end do 

57、    !至此，已经完成查找最大元素，查找完成以后与  第k行交换  !交换两行元素，其他不变    vtemp1=Ab(k,:)    vtemp2=Ab(id_max,:)           Ab(k,:)=vtemp2    Ab(id_max,:)=vtemp1   !以上一大段是为交换两行元素，交

58、换完成以后即按照消元法进行!#     do i=k+1,N       temp=Ab(i,k)/Ab(k,k)          Ab(i,:)=Ab(i,:)-temp*Ab(k,:)      end doend do!-! 经过上一步，Ab已经化为如下形式的矩阵!       &

59、#160;    | *  *  *  *  # |!     A b= | 0  *  *  *  # |!            | 0  0  *  *  # |!           

60、; | 0  0  0  *  # |!Aup(:,:)=Ab(1:N,1:N)bup(:)=Ab(:,N+1)!调用用上三角方程组的回带方法call uptri(Aup,bup,x,n)end subroutine solve subroutine uptri(A,b,x,N)!-subroutine  comment!  Version   :  V1.0    !  Coded by  :  syz !

61、0; Date      :  2010-4-8!-!  Purpose   :  上三角方程组的回带方法!                 Ax=b!-!  Input  parameters  :!       1.   A(

62、N,N)系数矩阵!       2.   b(N)右向量!       3.   N方程维数!  Output parameters  :!       1.  x  方程的根!       2.!  Common parameters  :!-implici

63、t real*8(a-z)integer:i,j,Nreal*8:A(N,N),b(N),x(N)x(N)=b(N)/A(N,N)!回带部分do i=n-1,1,-1       x(i)=b(i)   do j=i+1,N    x(i)=x(i)-a(i,j)*x(j)   end do    x(i)=x(i)/A(i,i)end doend subroutine uptriend module m_gauss prog

64、ram main!-program comment!  Version   :  V1.0    !  Coded by  :  syz !  Date      :  2010-4-8!-!  Purpose   :  高斯列主元消去法!    !-!  In put data  files :! 

65、;      1.  fin.txt  输入方程系数!       2.!  Output data files  :!       1. fout.txt  计算结果!       2.!-!  Post Script :!       1.

66、60;   需要准备输入数据!     !-use m_gaussimplicit real*8(a-z)integer,parameter: N=6integer:i,jreal*8:A(N,N),b(N),x(N) open(unit=11,file='fin.txt')open(unit=12,file='fout.txt')read(11,*)!读入A矩阵read(11,*)(A(i,j),j=1,N),i=1,N)!读入B向量read(11,*) bcall solve(A

67、,b,x,N)write(12,101)x101 format(T5,'高斯列主元消去法计算结果',/,T4,'x=',6(/F12.8)end program main6.追赶法module chase!-module coment!  Version     :  V1.0    !  Coded by    :  syz !  Date    

68、60;   :  2010-4-9!-!  Description :  三对角线方程组方法模块!    !-!  Post Script :!      1.!      2. !-contains subroutine solve(A,f,x,N)!-subroutine  comment!  Version   :  V1

69、.0    !  Coded by  :  syz !  Date      :  2010-4-9!-!  Purpose   :  追赶法计算三对角方程组!              Ax=f!-!  Input  parameters  :!&

70、#160;      1.  A系数矩阵!       2. f 右向量!  Output parameters  :!       1.  x方程的解!       2.  N维数!  Common parameters  :!-!  Post Script :!  

71、     1.   注意：该方法仅适用于三对角方程组!       2.!-  implicit real*8(a-z) integer:N real*8:A(N,N),f(N),x(N) real*8:L(2:N),u(N),d(1:N-1)  real*8:c(1:N-1),b(N),e(2:N) integer:i  real*8:y(N)!-把A矩阵复制给向量 e,b,

72、c do i=1,N    b(i)=a(i,i)end dodo i=1,N-1   c(i)=a(i,i+1)end dodo i=2,N  e(i)=a(i,i-1)end do!- do i=1,N-1 d(i)=c(i)end dou(1)=b(1)do i=2,N  L(i)=e(i)/u(i-1)  u(i)=b(i)-L(i)*c(i-1)end do !-开始回带,求得yy(1)=f(1)do i=2,N  y(i)=f(i)-L(i)*y(i-1

73、)end do!-开始回带，求得xx(n)=y(n)/u(n)do i=n-1,1,-1x(i)=(y(i)-c(i)*x(i+1)/u(i)end do end subroutine solveend module chaseprogram main!-program comment!  Version   :  V1.0    !  Coded by  :  syz !  Date      :

74、0; !-!  Purpose   :  追赶法计算三对角方程!    !-!  In put data  files :!       1.   fin.txt输入数据!       2.   fout.txt输出数据!  Output data files  :!   

75、60;   1.!       2.!-use chaseinteger,parameter:N=4real*8:A(N,N),f(N),x(N)open(unit=11,file='fin.txt')open(unit=12,file='fout.txt')read(11,*)read(11,*)(A(i,j),j=1,N),i=1,N)read(11,*)fcall solve(A,f,x,N)write(12,101)x101 format(T5,'追赶法计算结果',/,T4,