Part_II_數(shù)值分析程序設(shè)計.doc

數(shù)值分析程序設(shè)計Part IIFortran程序設(shè)計1 線性代數(shù)1.1 矩陣的加法和乘法矩陣作為數(shù)值分析的重要工具，在數(shù)值計算中具有非常重要的作用，在Fortran語言中利用數(shù)組來表示矩陣。在Fortran 90中可以直接對數(shù)組進(jìn)行計算，一個命令就可以完成矩陣的加法、減法。Real a(m,n), b(m,n),c(m,n)C=a+b ！矩陣加法C=a-b ！矩陣減法在Fortran 77中必須利用循環(huán)完成矩陣的加減法：do r=1,m do c=1,n c(r,c)=a(r,c)+b(r,c) end doend do矩陣的乘法，在Fortran 90中不能直接用乘號來做，必須調(diào)用庫函數(shù)。C=a*b ！這個命令是做c(I,j)=a(I,j)*b(I,j)C=matmul(a,b) ！庫函數(shù)matmul可以做矩陣相乘在Fortran 77中則要自己寫矩陣乘法程序代碼進(jìn)行計算，下面是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 用字符串來設(shè)定輸出格式 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 ! 循環(huán)的計數(shù)器 ! 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，那么一定要避免下面的情況。使用下面的方法，在子程序中會讀到錯誤的矩陣內(nèi)容。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實際聲明數(shù)組大小為10*10，傳遞到子程序中卻把它聲明為2*2，會有下面的對應(yīng)情況發(fā)生：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)這是由于多維數(shù)組的存儲方式?jīng)Q定的，在數(shù)據(jù)讀取時應(yīng)當(dāng)注意。1.2 三角矩陣通過矩陣的初等變換，可以將矩陣化成上三角矩陣、下三角矩陣或者對角矩陣。下面程序?qū)⒕仃嚪謩e化成上三角矩陣和下三角矩陣：=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) ! 用字符串來設(shè)定輸出格式 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的功能可以少一層循環(huán) 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的功能可以少一層循環(huán) 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程序執(zhí)行結(jié)果：矩陣的三角化可以應(yīng)用于行列式計算、求解逆矩陣、線性方程組求解等方面。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的功能可以少一層循環(huán) 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消去法是直接求解線性方程組的數(shù)值方法。通過矩陣的初等變換，將矩陣化成對角矩陣求解。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,及數(shù)組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,可以中止斷行發(fā)生,使下一次的 ! write接續(xù)在同一行當(dāng)中. 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! 求解聯(lián)立式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程序執(zhí)行結(jié)果：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,及數(shù)組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,可以中止斷行發(fā)生,使下一次的 ! write接續(xù)在同一行當(dāng)中. 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) ! 調(diào)用排序的子程序 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之后的數(shù)字不用比較 ! 如果A(J) A(J+1) 就把這兩個數(shù)值交換 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個數(shù)據(jù)中最小的一個，把它和數(shù)列的第一個數(shù)字交換位置；（2）找出剩余n-1個數(shù)據(jù)中最小的一個，把它和數(shù)列的第二個數(shù)字交換位置；（3）找出剩余n-2個數(shù)據(jù)中最小的一個，把它和數(shù)列的第三個數(shù)字交換位置；（4）（5）一直做到只剩一個數(shù)據(jù)為止。計算程序如下：! 選擇排序法范例! 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/) ! 排序的數(shù)據(jù) write(*,(Source=,10I3) A call SELECTION_SORT(A,N) ! 調(diào)用排序的子程序 write(*,(Sort=,10I3) A stopend program! 選擇排序法的子程序subroutine SELECTION_SORT(A,N) implicit none integer : N,A(N) integer I,J ! 循環(huán)計數(shù)器 integer MIN ! 找出每一輪中的最小值 integer TEMP ! 交換數(shù)據(jù)時使用 do I=1,N MIN=A(I) ! 暫時令A(yù)(I)是最小值 do J=I+1,N if ( MIN A(J) ) then ! 發(fā)現(xiàn)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 順序搜索法順序搜索是一種最簡單的搜索方法，其原理是將數(shù)據(jù)一個一個拿出來，看他是不是我們所需要的數(shù)據(jù)。順序搜索法的程序如下：! 順序查找法范例! 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/) ! 存放數(shù)據(jù)組的類型 integer KEY ! 記錄所要找的值 integer LOC integer, external : SEQUENTIAL_SEARCH write(*,(Source=,10I3) A write(*,*) Input KEY: read (*,*) KEY ! 鍵入待尋數(shù)據(jù) ! 調(diào)用順序查找的函數(shù) 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 ! 循環(huán)的計數(shù)器 do I=1,N ! 開始做掃瞄, 最多做N次. if ( KEY=A(I) ) then ! 找到了, 返回數(shù)字在類型中的位置 SEQUENTIAL_SEARCH=I returnend if end do ! 沒找到時返回-1 SEQUENTIAL_SEARCH=0 returnend function4.2.2 二元搜索二元搜索法必須配合排序好的數(shù)據(jù)才能使用，假設(shè)所要尋找的數(shù)據(jù)值為K，數(shù)據(jù)存放在數(shù)組中，搜索的步驟為：（1）取出數(shù)組的中間值M和K來比較，如果M=K，結(jié)束。如果KM，因為數(shù)組早已作好排序，所以K值一定是在數(shù)組的后半段。如果K,10I3) A write(*,*) Input KEY: read (*,*) KEY ! 調(diào)用查找的子程序 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 ! 記錄每一個小組的類型結(jié)束位置 integer M ! 記錄每一個小組的類型中間位置 ! 一開始的小組范圍就是整個類型 L=1 R=