Part_II_数值分析程序设计.doc

数值分析程序设计Part IIFortran程序设计1 线性代数1.1 矩阵的加法和乘法矩阵作为数值分析的重要工具，在数值计算中具有非常重要的作用，在Fortran语言中利用数组来表示矩阵。在Fortran 90中可以直接对数组进行计算，一个命令就可以完成矩阵的加法、减法。Real a(m,n), b(m,n),c(m,n)C=a+b ！矩阵加法C=a-b ！矩阵减法在Fortran 77中必须利用循环完成矩阵的加减法：do r=1,m do c=1,n c(r,c)=a(r,c)+b(r,c) end doend do矩阵的乘法，在Fortran 90中不能直接用乘号来做，必须调用库函数。C=a*b ！这个命令是做c(I,j)=a(I,j)*b(I,j)C=matmul(a,b) ！库函数matmul可以做矩阵相乘在Fortran 77中则要自己写矩阵乘法程序代码进行计算，下面是Fortran 77固定格式编写的矩阵相乘程序代码：CC 矩阵乘法范例C By Perng 1997/9/17 PROGRAM MATMUL_DEMO IMPLICIT NONE INTEGER N PARAMETER(N=3) INTEGER A(N,N) ! MATRIX A INTEGER B(N,N) ! MATRIX B INTEGER C(N,N) ! MATRIX C DATA B /1,2,3,4,5,6,7,8,9/ DATA C /9,8,7,6,5,4,3,2,1/ CALL MATMUL(A,B,N,N,C,N,N) WRITE(*,*) Matrix A: CALL OUTPUT(A,N) STOP ENDCC 输出矩阵的子程序C SUBROUTINE OUTPUT(A,N) IMPLICIT NONE INTEGER N,A(N,N) INTEGER I,J CHARACTER FOR*20 DATA FOR /(?(1X,I3)/C 用字符串来设定输出格式 WRITE( FOR(2:3), (I2) ) N DO I=1,N WRITE( *, FMT=FOR ) (A(I,J),J=1,N) END DO RETURN ENDCC 矩阵乘法的子程序C SUBROUTINE MATMUL(A,B,BR,BC,C,CR,CC) IMPLICIT NONE INTEGER BR ! Row of Matrix B INTEGER BC ! Column of Matrix B INTEGER B(BR,BC) ! Matrix B INTEGER CR ! Row of Matrix C INTEGER CC ! Column of Matrix C INTEGER C(CR,CC) ! Matrix C INTEGER A(BR,CC) ! Matrix A INTEGER I,J,K ! 循环的计数器 ! BC若不等于CR, 这两个矩阵无法相乘 IF ( BC .NE. CR ) THEN WRITE(*,*) Matrix size error! STOP END IF DO I=1,BR DO J=1,CC A(I,J)=0 DO K=1,BC A(I,J)=A(I,J)+B(I,K)*C(K,J) END DO END DO END DO RETURN END需要注意的是：假如矩阵一开始声明了10*10的大小，而所需要使用的矩阵大小只有2*2，那么一定要避免下面的情况。使用下面的方法，在子程序中会读到错误的矩阵内容。Program mianImplicit noneInteger A(10,10)A(1,1)=1.0A(2,1)=2.0A(1,2)=1.0A(2,2)=2.0Call sub(A,2,2).end program mainsubroutine sub(matrix, row, col)implicit noneinteger row, colinteger matrix(row, col)end subroutine sub实际声明数组大小为10*10，传递到子程序中却把它声明为2*2，会有下面的对应情况发生：matrix(1,1)=a(1,1)=1.0matrix(2,1)=a(2,1)=2.0matrix(1,2)=a(3,1)=0.0/=a(1,2)matrix(2,2)=a(4,1)=0.0/=a(2,2)这是由于多维数组的存储方式决定的，在数据读取时应当注意。1.2 三角矩阵通过矩阵的初等变换，可以将矩阵化成上三角矩阵、下三角矩阵或者对角矩阵。下面程序将矩阵分别化成上三角矩阵和下三角矩阵：=UPPERLOWER.F90module LinearAlgebra implicit nonecontains! 输出矩阵的子程序subroutine output(matrix) implicit none integer : m,n real : matrix(:,:) integer : i character(len=20) : for=(?(1x,f6.3) m = size(matrix,1) n = size(matrix,2) ! 用字符串来设定输出格式 write( FOR(2:3), (I2) ) N do i=1,Nwrite( *, FMT=FOR ) matrix(i,:) end do returnend subroutine output! 求上三角矩阵的子程序subroutine Upper(matrix) implicit none real : matrix(:,:) integer : M,N integer : I,J real : E M=size(matrix,1) N=size(matrix,2) do I=1,N-1do J=I+1,M E=matrix(J,I)/matrix(I,I) ! 用90的功能可以少一层循环 matrix(J,I:M)=matrix(J,I:M)-matrix(I,I:M)*Eend do end do returnend subroutine Upper! 求下三角矩阵的子程序subroutine Lower(matrix) implicit none real : matrix(:,:) integer : M,N real : I,J,E M = size(matrix,1) N = size(matrix,2) do I=N,2,-1 do J=I-1,1,-1 E=matrix(J,I)/matrix(I,I) ! 用90的功能可以少一层循环 matrix(J,1:I)=matrix(J,1:I)-matrix(I,1:I)*E end do end do returnend subroutine Lowerend moduleprogram main use LinearAlgebra implicit none integer, parameter : N = 3! Size of Matrix real : A(N,N) = reshape( (/1,2,1,3,2,3,2,3,4/),(/N,N/) ) real : B(N,N) write(*,*) Matrix A: call output(A) B=A write(*,*) Upper: call Upper(B) call output(B) B=A write(*,*) Lower: call Lower(B) call output(B) stopend program程序执行结果：矩阵的三角化可以应用于行列式计算、求解逆矩阵、线性方程组求解等方面。1.3 矩阵的行列式计算方阵的行列式计算，采用矩阵三角化方法，矩阵的行列式等于其三角化矩阵对角线元素的乘积。下面程序先把矩阵化成上三角矩阵，再求行列式的值。DETERMINANT.F90module LinearAlgebra implicit nonecontains! 求矩阵的Determinant值real function Determinant(matrix) real : matrix(:,:) real, allocatable : ma(:,:) integer : i,N N = size(matrix,1) allocate(ma(N,N) ma = matrix call Upper(ma) Determinant = 1.0 do i=1,N Determinant = Determinant*ma(i,i) end doend function! 求上三角矩阵的子程序subroutine Upper(matrix) real : matrix(:,:) integer : M,N integer : I,J real : E M=size(matrix,1) N=size(matrix,2) do I=1,N-1do J=I+1,M E=matrix(J,I)/matrix(I,I) ! 用90的功能可以少一层循环 matrix(J,I:M)=matrix(J,I:M)-matrix(I,I:M)*Eend do end do returnend subroutine Upperend moduleprogram main use LinearAlgebra implicit none integer, parameter : N = 3! Size of Matrix real : A(N,N) = reshape( (/1,2,1,3,2,3,2,3,4/),(/N,N/) ) write(*,(det(A)=,F6.2) Determinant(A) stopend program1.4 Gauss-Jordan消去法求解线性方程组线性方程组的求解可以分为直接法和迭代法。直接法有Gauss消去法和直接三角分解法。迭代法有Jacobi迭代法和Gauss-Seidel迭代法等。Gauss-Jordan消去法是直接求解线性方程组的数值方法。通过矩阵的初等变换，将矩阵化成对角矩阵求解。GAUSS-JORDAN.F90module LinearAlgebra implicit nonecontains! Gauss_Jordan法subroutine Gauss_Jordan(A,S,ANS) implicit none real : A(:,:) real : S(:) real : ANS(:) real, allocatable : B(:,:) integer : i, N N = size(A,1) allocate(B(N,N) ! 保存原先的矩阵A,及数组S B=A ANS=S ! 把B化成对角线矩阵(除了对角线外,都为0) call Upper(B,ANS,N) ! 先把B化成上三角矩阵 call Lower(B,ANS,N) ! 再把B化成下三角矩阵 ! 求解 forall(i=1:N) ANS(i)=ANS(i)/B(i,i) end forall returnend subroutine Gauss_Jordan! 输出等式subroutine output(M,S) implicit none real : M(:,:), S(:) integer : N,i,j N = size(M,1) ! write中加上advance=no,可以中止断行发生,使下一次的 ! write接续在同一行当中. do i=1,N write(*,(1x,f5.2,a1), advance=NO) M(i,1),A do j=2,N if ( M(i,j) 0 ) then write(*,(-,f5.2,a1),advance=NO) -M(i,j),char(64+j) else write(*,(+,f5.2,a1),advance=NO) M(i,j),char(64+j) end if end do write(*,(=,f8.4) S(i) end do returnend subroutine output! 求上三角矩阵的子程序subroutine Upper(M,S,N) implicit none integer : N real : M(N,N) real : S(N) integer : I,J real : E do I=1,N-1 do J=I+1,N E=M(J,I)/M(I,I) M(J,I:N)=M(J,I:N)-M(I,I:N)*E S(J)=S(J)-S(I)*E end do end do returnend subroutine Upper! 求下三角矩阵的子程序subroutine Lower(M,S,N) implicit none integer : N real : M(N,N) real : S(N) integer : I,J real : E do I=N,2,-1 do J=I-1,1,-1 E=M(J,I)/M(I,I) M(J,1:N)=M(J,1:N)-M(I,1:N)*E S(J)=S(J)-S(I)*E end do end do returnend subroutine Lowerend module! 求解联立式program main use LinearAlgebra implicit none integer, parameter : N=3 ! Size of Matrix real : A(N,N)=reshape( (/1,2,3,4,5,6,7,8,8/),(/N,N/) ) real : S(N)=(/12,15,17/) real : ans(N) integer : i write(*,*) Equation: call output(A,S) call Gauss_Jordan(A,S,ANS) write(*,*) Ans: do i=1,N write(*,(1x,a1,=,F8.4) char(64+i),ANS(i) end do stopend program程序执行结果：1.5 逆矩阵求解逆矩阵的算法与Gauss-Jordan消去法求解方程类似，在矩阵的右边添加一个同阶的单位矩阵，将左边矩阵变换为单位矩阵，右边的单位矩阵变换为矩阵的逆矩阵。程序如下：INVERSE.F90module LinearAlgebra implicit nonecontains! Gauss_Jordan法subroutine Gauss_Jordan(A,S,ANS) implicit none real : A(:,:) real : S(:) real : ANS(:) real, allocatable : B(:,:) integer : i, N N = size(A,1) allocate(B(N,N) ! 保存原先的矩阵A,及数组S B=A ANS=S ! 把B化成对角线矩阵(除了对角线外,都为0) call Upper(B,ANS,N) ! 先把B化成上三角矩阵 call Lower(B,ANS,N) ! 再把B化成下三角矩阵 ! 求解 forall(i=1:N) ANS(i)=ANS(i)/B(i,i) end forall returnend subroutine Gauss_Jordan! 输出等式subroutine output(M,S) implicit none real : M(:,:), S(:) integer : N,i,j N = size(M,1) ! write中加上advance=no,可以中止断行发生,使下一次的 ! write接续在同一行当中. do i=1,N write(*,(1x,f5.2,a1), advance=NO) M(i,1),A do j=2,N if ( M(i,j) ,10I3) A call BUBBLE_SORT(A,N) ! 调用排序的子程序 write(*,(Sort=,10I3) A stopend programsubroutine BUBBLE_SORT(A,N) implicit none integer : N,A(N) integer I,J, TEMP do I=N-1,1,-1 ! 开始做N-1次的扫瞄 do J=1,I ! 一对一对的来比较，I之后的数字不用比较 ! 如果A(J) A(J+1) 就把这两个数值交换 if ( A(J) A(J+1) ) then TEMP=A(J) A(J)=A(J+1) A(J+1)=TEMP end if end do end do returnend subroutine4.1.2 选择排序法选择排序法的基本原理为：（1）找出全部n个数据中最小的一个，把它和数列的第一个数字交换位置；（2）找出剩余n-1个数据中最小的一个，把它和数列的第二个数字交换位置；（3）找出剩余n-2个数据中最小的一个，把它和数列的第三个数字交换位置；（4）（5）一直做到只剩一个数据为止。计算程序如下：! 选择排序法范例! By Perng 1997/8/29program SELECTION_SORT_DEMO implicit none integer, parameter : N=10 integer : A(N)=(/6,2,8,4,0,9,3,5,1,7/) ! 排序的数据 write(*,(Source=,10I3) A call SELECTION_SORT(A,N) ! 调用排序的子程序 write(*,(Sort=,10I3) A stopend program! 选择排序法的子程序subroutine SELECTION_SORT(A,N) implicit none integer : N,A(N) integer I,J ! 循环计数器 integer MIN ! 找出每一轮中的最小值 integer TEMP ! 交换数据时使用 do I=1,N MIN=A(I) ! 暂时令A(I)是最小值 do J=I+1,N if ( MIN A(J) ) then ! 发现A(I)不是最小 TEMP=A(J) ! 把A(I)、A(J)交换 A(J)=A(I) A(I)=TEMP MIN=A(I) end ifend do end do returnend subroutine 4.2 搜索4.2.1 顺序搜索法顺序搜索是一种最简单的搜索方法，其原理是将数据一个一个拿出来，看他是不是我们所需要的数据。顺序搜索法的程序如下：! 顺序查找法范例! By Perng 1997/8/31program SEQUENTIAL_SEARCH_DEMO implicit none integer, parameter : N=10 integer : A(N)=(/6,2,8,4,0,9,3,5,1,7/) ! 存放数据组的类型 integer KEY ! 记录所要找的值 integer LOC integer, external : SEQUENTIAL_SEARCH write(*,(Source=,10I3) A write(*,*) Input KEY: read (*,*) KEY ! 键入待寻数据 ! 调用顺序查找的函数 LOC = SEQUENTIAL_SEARCH(A,N,KEY) if ( LOC/=0 ) then write(*,(A(,I2, )=I3) LOC,KEY else write(*,*) Not found end if stopend program! 顺序查找法的子程序integer function SEQUENTIAL_SEARCH(A,N,KEY) implicit none integer N, A(N) integer KEY ! 所要寻找的值 integer I ! 循环的计数器 do I=1,N ! 开始做扫瞄, 最多做N次. if ( KEY=A(I) ) then ! 找到了, 返回数字在类型中的位置 SEQUENTIAL_SEARCH=I returnend if end do ! 没找到时返回-1 SEQUENTIAL_SEARCH=0 returnend function4.2.2 二元搜索二元搜索法必须配合排序好的数据才能使用，假设所要寻找的数据值为K，数据存放在数组中，搜索的步骤为：（1）取出数组的中间值M和K来比较，如果M=K，结束。如果KM，因为数组早已作好排序，所以K值一定是在数组的后半段。如果K,10I3) A write(*,*) Input KEY: read (*,*) KEY ! 调用查找的子程序 LOC=BINARY_SEARCH(A,N,KEY) if ( LOC/=0 ) then write(*,(A(,I2, )=I3) LOC,KEY else write(*,*) Not found end if stopend program! 折半查找法的子程序integer function BINARY_SEARCH(A,N,KEY) implicit none integer N,A(N) integer KEY ! 所要寻找的值 integer L ! 记录每一个小组的类型起始位置 integer R ! 记录每一个小组的类型结束位置 integer M ! 记录每一个小组的类型中间位置 ! 一开始的小组范围就是整个类型 L=1 R=