《计算方法》Fortran版

1、第一章：线性方程组的解法一 、 矩阵分解与线性方程组直接方法1. LU 分解解方程module LU!-module coment!Version:V1.0!Coded by:syz!Date:!-! Description : LU 分解解方程!-containssubroutine 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)calldowntri(L,b,y,N)call uptri(U,y,x

2、,N)end subroutine solvesubroutine doolittle(A,L,U,N)!-subroutinecomment!Version:V1.0!Coded by:syz!Date:!-!Purpose:LU 分解之 Doolittle函数!A=LU!-!Input parameters :1/54! 1. A方阵! 2. N阶数!Output parameters:! 1. L! 2. U!Common parameters:!-! Post Script :! 1.! 2.!-implicit real*8(a-z)integer:N,i,k,rreal*8:A(N

3、,N),L(N,N),U(N,N)!U 的第一行U(1,:)=A(1,:)!L 的第一列L(:,1)=a(:,1)/U(1,1)do k=2,Nl(k,k)=1do j=k,ns=0do m=1,k-1s=s+l(k,m)*u(m,j)end dou(k,j)=a(k,j)-send dodo i=k+1,ns=0do m=1,k-1s=s+l(i,m)*u(m,k)end dol(i,k)=(a(i,k)-s)/u(k,k)end do2/54end doend subroutine doolittlesubroutine uptri(A,b,x,N)!-subroutinecomment!

4、Version:V1.0!Coded by:syz!Date:2010-4-8!-!Purpose:上三角方程组的回带方法!Ax=b!-!Input parameters :! 1. A(N,N) 系数矩阵! 2. b(N) 右向量! 3. N 方程维数!Output parameters:! 1. x 方程的根! 2.!Common parameters:!-implicit real*8(a-z)integer:i,j,k,Nreal*8:A(N,N),b(N),x(N)x(N)=b(N)/A(N,N)! 回带部分do i=n-1,1,-1x(i)=b(i)do j=i+1,Nx(i)=x

5、(i)-a(i,j)*x(j)end dox(i)=x(i)/A(i,i)end doend subroutine uptri3/54subroutine downtri(A,b,x,N)!-subroutinecomment!Version:V1.0!Coded by:syz!Date:2010-4-9!-!Purpose:下三角方程组的回带方法!Ax=b!-!Input parameters :! 1. A(N,N) 系数矩阵! 2. b(N) 右向量! 3. N 方程维数!Output parameters:! 1. x 方程的根! 2.!Common parameters:!-impl

6、icit 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,Nx(k)=b(k)do i=1,k-1x(k)=x(k)-a(k,i)*x(i)end dox(k)=x(k)/a(k,k)end doend subroutine downtriend module LUprogram main!-program comment!Version:V1.0!Coded by:syz!Date:2010-4-84/54!-!Purpose:LU 分解计算线性方程组!Ax=b!-!In put datafiles

7、:!1.A,b! 2.!Output data files:!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,N)write(12,101)x101 format(T5,'LU分解计算线性方程组计算

8、结果',/,4(/,F10.6)end program main2. LU 分解之 Crout module crout containssubroutine solve(A,L,U,N)!-subroutinecomment!Version:V1.0!Coded by:syz!Date:!-!Purpose:LU 之 Crout分解!A=LU!-5/54!Inputparameters:! 1. A方阵! 2. N阶数!Output parameters:! 1. L! 2. U!Common parameters:!-! Post Script :! 1.! 2.!-implici

9、t 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,Ndo i=k,ns=0do r=1,k-1s=s+l(i,r)*u(r,k)end dol(i,k)=a(i,k)-send dodo j=k+1,ns=0do r=1,k-1s=s+l(k,r)*u(r,j)end dou(k,j)=(a(k,j)-s)/l(k,k)end dou(k,k)=16/54end doend subroutine solveend module cr

10、outprogram main!-program comment!Version:V1.0!Coded by :syz!Date:2010-4-8!-!Purpose:Crot 分解!-!In put datafiles :!1.A,N! 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(11,*)read(11,

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,Nwrite(12,22)L(i,:)end do22 format(4F10.6)7/54! 输出U矩阵write(12,*)'U='do i=1,Nwrite(12,22)U(i,:)end do23 format(4F10.6) end program main3. LU 分解之 Doolittle module d

12、oolittle!-module coment!Version:V1.0!Coded by: syz!Date:!-!Description :LU 分解之 doolittle模块!-containssubroutine solve(A,L,U,N)!-subroutinecomment!Version:V1.0!Coded by:syz!Date:!-!Purpose:LU 分解之 Doolittle函数!A=LU!-!Input parameters :! 1. A方阵! 2. N阶数!Output parameters:! 1. L! 2. U!Common parameters:!-!

13、 Post Script :! 1.! 2.!-8/54implicit 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)/U(1,1)do k=2,Nl(k,k)=1do j=k,ns=0do m=1,k-1s=s+l(k,m)*u(m,j)end dou(k,j)=a(k,j)-send dodo i=k+1,ns=0do m=1,k-1s=s+l(i,m)*u(m,k)end dol(i,k)=(a(i,k)-s)/u(k,k)end doend do

14、end subroutine solveend module doolittleprogram main!-program comment!Version:V1.09/54!Coded by :syz!Date:2010-4-8!-!Purpose:Doolittle 分解!-!In put datafiles :!1.A,N! 2.!Output 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(

15、unit=12,file='fout.txt')read(11,*)read(11,*)(A(i,j),j=1,N),i=1,N)call solve(A,L,U,N)write(12,21)21 format(T10,'Doolittle分解 ',/)! 输出 L矩阵write(12,*)'L='do i=1,Nwrite(12,22)L(i,:)end do22 format(3F10.6)! 输出U矩阵write(12,*)'U='do i=1,Nwrite(12,22)U(i,:)end do23 format(3F10.

16、6) end program main4. 高斯消去法10/54module gauss!-module coment!Version:V1.0!Coded by:syz!Date:2010-4-8!-! Description : 高斯消去法模块!-!Contains:!1.solve方法函数!2.!-containssubroutine solve(A,b,x,N)!-subroutinecomment!Version:V1.0!Coded by:syz!Date:2010-4-8!-!Purpose:高斯消去法!Ax=b!-!Input parameters :! 1. A(N,N) 系

17、数矩阵! 2. 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:Aup(N,N),bup(N)!Ab 为增广矩阵Abreal*8:Ab(N,N+1)11/54Ab(1:N,1:N)=AAb(:,N+1)=b!-! 这段是 高斯消去法的核心部分do k=1,N-1do i=k+1,Ntemp=Ab(i,k)/Ab(k,k)Ab(i,:)=Ab(i,:)-temp*Ab(k,

18、:)end doend do!-! 经过上一步， Ab 已经化为如下形式的矩阵!| *# |!A b= | 0* # |!| 00*# |!| 000*# |!Aup(:,:)=Ab(1:N,1:N)bup(:)=Ab(:,N+1)! 调用用上三角方程组的回带方法call uptri(Aup,bup,x,n)end subroutine solvesubroutine uptri(A,b,x,N)!-subroutinecomment!Version:V1.0!Coded by:syz!Date:2010-4-8!-!Purpose:上三角方程组的回带方法!Ax=b!-!Inputparame

19、ters :!1.A(N,N) 系数矩阵12/54! 2. b(N) 右向量! 3. N 方程维数!Output parameters:! 1. x 方程的根! 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,-1x(i)=b(i)do j=i+1,Nx(i)=x(i)-a(i,j)*x(j)end dox(i)=x(i)/A(i,i)end doend subroutine uptriend module gaussp

20、rogram main!-program comment!Version:V1.0!Coded by:syz!Date:2010-4-8!-!Purpose:高斯消去法!-!In put datafiles :! 1. fin.txt 输入方程系数! 2.!Output data files:!1. fout.txt计算结果! 2.13/54!-! 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=1

21、1,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,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

22、:syz!Date:2010-4-8!-! Description : 高斯列主元消去法模块!-!Contains:!1.solve方法函数!2.!-14/54containssubroutine solve(A,b,x,N)!-subroutinecomment!Version:V1.0!Coded by:syz!Date:2010-4-8!-!Purpose:高斯列主元消去法!Ax=b!-!Input parameters :! 1. A(N,N) 系数矩阵! 2. b(N) 右向量! 3. N 方程维数!Output parameters:! 1. x 方程的根! 2.!Common p

23、arameters:!-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)real*8:vtemp1(N+1),vtemp2(N+1)Ab(1:N,1:N)=AAb(:,N+1)=b!# #! 这段是 列主元消去法的核心部分do k=1,N-1elmax=dabs(Ab(k,k)id_max=k15/54!这段为查找主元素!这段程序的主要目的不是为了赋值最大元素给elmax ，而是为了找出最大元

24、素对应的标号do i=k+1,nif (dabs(Ab(i,k)>elmax) thenelmax=Ab(i,k)id_max=iend ifend do! 至此，已经完成查找最大元素，查找完成以后与第 k 行交换! 交换两行元素，其他不变 vtemp1=Ab(k,:) vtemp2=Ab(id_max,:)Ab(k,:)=vtemp2Ab(id_max,:)=vtemp1! 以上一大段是为交换两行元素，交换完成以后即按照消元法进行!# #do i=k+1,Ntemp=Ab(i,k)/Ab(k,k)Ab(i,:)=Ab(i,:)-temp*Ab(k,:)end doend do!-! 经

25、过上一步， Ab 已经化为如下形式的矩阵!| *# |!A b= | 0* # |!| 00*# |!| 000*# |!Aup(:,:)=Ab(1:N,1:N)bup(:)=Ab(:,N+1)16/54! 调用用上三角方程组的回带方法call uptri(Aup,bup,x,n)end subroutine solvesubroutine uptri(A,b,x,N)!-subroutinecomment!Version:V1.0!Coded by:syz!Date:2010-4-8!-!Purpose:上三角方程组的回带方法!Ax=b!-!Input parameters :! 1. A(

26、N,N) 系数矩阵! 2. b(N) 右向量! 3. N 方程维数!Output parameters:! 1. x 方程的根! 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,-1x(i)=b(i)do j=i+1,Nx(i)=x(i)-a(i,j)*x(j)end dox(i)=x(i)/A(i,i)end doend subroutine uptriend module m_gauss17/54program ma

27、in!-program comment!Version:V1.0!Coded by:syz!Date:2010-4-8!-!Purpose:高斯列主元消去法!-!In put datafiles :! 1. fin.txt 输入方程系数! 2.!Output data files:!1. fout.txt计算结果!2.!-!Post Script :!1.需要准备输入数据!-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='f

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

29、2010-4-9!-! Description : 三对角线方程组方法模块!-! Post Script :! 1.! 2.!-containssubroutine solve(A,f,x,N)!-subroutinecomment!Version:V1.0!Coded by:syz!Date:2010-4-9!-!Purpose:追赶法计算三对角方程组!Ax=f!-!Input parameters :! 1. A 系数矩阵! 2. f 右向量!Output parameters:! 1. x 方程的解! 2. N维数!Common parameters:!-! Post Script :! 1. 注意：该方法仅适用于三对角方程组! 2.!-implicit real*8(a-z)integer:N19/54real*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:ireal*8:y(N)!-把 A 矩阵复制给向量e,b,cdo i=1,Nb(i)=a(i,i)end dodo i=1,N-1c(i)=a(i,i+1)end dodo i=2,Ne(i)=a(i,i-1)end do!-do i=1,N-1d(i)=c(i)end dou(1