计算方法Fortran版.docx

1、计算方法Fortran版-CAL-FENGHAI.-(YICAI)-Company One1 -CAL-本页仅作为文档封面，使用请直接删除计算方法Fortran版(总20页)第一章：线性方程组的解法一 、矩阵分解与线性方程组直接方法1. LU分解解方程module LU!-module coment! Version : ! Coded by : syz! Date :!-! Description : LU 分解解方程!-contains?subroutine solve(A,b,x,N)implicit real*8(a-z)integer:Nreal*8:A(N,N),b(N),x(N)

2、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 doolittle(A,L,U,N)!-subroutine comment! Version : ! Coded by : syz! Date :!-! Purpose : LU分解之Doolittle函数! A=LU!-! Input parameters :! 1. A 方阵! 2. N 阶数! Output parameters :! 1.

3、 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)/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 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

4、 do l(i,k)=(a(i,k)-s)/u(k,k) end doend doend subroutine doolittlesubroutine uptri(A,b,x,N)!-subroutine comment! Version : ! 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

5、 :!-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,-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 comment! Version : ! Coded by : syz! Date : 2010-4-9!-! Purpose : 下三角方程

6、组的回带方法! 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,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) end do x(k)=x(k)/a(k,k)end doend subroutine do

7、wntriend module LU?program main!-program comment! Version : ! Coded by : syz! Date : 2010-4-8!-! Purpose : LU分解计算线性方程组! Ax=b!-! In put data files :! 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=)open(unit=12,file=)

8、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分解计算线性方程组计算结果,LU分解之Croutmodule croutcontains?subroutine solve(A,L,U,N)!-subroutine comment! Version : ! Coded by : syz! Date :!-! Purpose : LU之Crout分解! A=LU!-! Input parameters :! 1. A 方阵! 2. N 阶数! O

9、utput parameters :! 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)!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(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

10、,r)*u(r,j) end do u(k,j)=(a(k,j)-s)/l(k,k) end do u(k,k)=1end doend subroutine solveend module crout?program main!-program comment! Version : ! Coded by : syz! Date : 2010-4-8!-! Purpose : Crot 分解!-! In put data files :! 1. A,N! 2.! Output data files :! 1. L,U! 2.!-use croutinteger,parameter:N=4real

11、*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)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!输出U矩阵write(12,*)U=do i=1,N?write(12,22)U(i,:)end do23 formatend program main3.

12、 LU分解之Doolittlemodule doolittle!-module coment! Version : ! Coded by : syz! Date :!-! Description : LU 分解之doolittle模块!-contains?subroutine solve(A,L,U,N)!-subroutine comment! Version : ! Coded by : syz! Date :!-! Purpose : LU分解之Doolittle函数! A=LU!-! Input parameters :! 1. A 方阵! 2. N 阶数! Output parame

13、ters :! 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,:)=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) 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)

14、*u(m,k) end do l(i,k)=(a(i,k)-s)/u(k,k) end doend doend subroutine solveend module doolittle?program main!-program comment! Version : ! Coded by : syz! Date : 2010-4-8!-! Purpose : Doolittle 分解!-! In put data files :! 1. A,N! 2.! Output data files :! 1. L,U! 2.!-use doolittleinteger,parameter:N=3rea

15、l*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)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!输出U矩阵write(12,*)U=do i=1,N?write(12,22)U(i,:)end do23 formatend program main

16、4.高斯消去法module gauss!-module coment! Version : ! Coded by : syz! Date : 2010-4-8!-! Description : 高斯消去法模块!-! Contains :! 1. solve 方法函数! 2.!-containssubroutine solve(A,b,x,N)!-subroutine comment! Version : ! 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)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,:) e

18、nd doend do!-! 经过上一步，Ab已经化为如下形式的矩阵! | * * * * # |! A b= | 0 * * * # |! | 0 0 * * # |! | 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 : ! Coded by : syz! Date : 2010-4-8!-! Purpose : 上

19、三角方程组的回带方法! 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,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

20、 subroutine uptriend module gaussprogram main!-program comment! Version : ! Coded by : syz! Date : 2010-4-8!-! Purpose : 高斯消去法!-! In put data files :! 1. 输入方程系数! 2.! Output data files :! 1. 计算结果! 2.!-! Post Script :! 1. 需要准备输入数据!-use gaussimplicit real*8(a-z)integer,parameter: N=4integer:i,jreal*8:A

21、(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,*) bcall solve(A,b,x,N)write(12,101)x101 format(T5,高斯消去法计算结果,/,T4,x=,4(/)end program main5.列主元消去法module m_gauss!-module coment! Version : ! Coded by : syz! Date : 2010-4-8!-! Description : 高斯列主元消去法模块!-! Contains :! 1. solve 方法函数! 2.!-containssubroutine solve(A,b,x,N)!-subroutine comment! Version : ! Coded by : syz! Date : 2010-4-8!-! Purpose : 高斯列主元消去法! Ax=b!-! Input p