PartII数值分析程序设计

1、数值分析程序设计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

2、=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、,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)' )

4、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) ! Matr

5、ix 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，那么一定要避

6、免下面的情况。使用下面的方法，在子程序中会读到错误的矩阵内容。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,

7、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

8、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) impli

9、cit 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

10、 : 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, para

11、meter : 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程序执行结果：矩阵的

12、三角化可以应用于行列式计算、求解逆矩阵、线性方程组求解等方面。1.3 矩阵的行列式计算方阵的行列式计算，采用矩阵三角化方法，矩阵的行列式等于其三角化矩阵对角线元素的乘积。下面程序先把矩阵化成上三角矩阵，再求行列式的值。DETERMINANT.F90module LinearAlgebra implicit nonecontains! 求矩阵的Determinant值real function Determinant(matrix) real : matrix(:,:) real, allocatable : ma(:,:) integer : i,N N = size(matrix,1) al

13、locate(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

14、的功能可以少一层循环 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)") Deter

15、minant(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

16、 : 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 Gaus

17、s_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(*,"

18、;('-',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 n

19、one 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

20、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)=(

21、/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消去法求解方程类似，在矩阵的右边添加一个同阶的单位矩阵，将左边

22、矩阵变换为单位矩阵，右边的单位矩阵变换为矩阵的逆矩阵。程序如下：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化成对角线矩阵(除了对

23、角线外,都为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",可以中止断行发生,使下一次的

24、! 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)

25、,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

26、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! 求解联立式prog

27、ram 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 w

28、rite(*,"(1x,a1,'=',F8.4)") char(64+i),ANS(i) end do stopend program2 数值积分对于积分进行数值积分，就是利用被积函数在区间上的一些节点处的函数值的线性组合去近似。其中称为积分节点，称为积分系数，为近似求积公式。积分节点和积分系数的不同决定了不同的积分公式。常用的有梯形积分公式、Simpson积分公式和Gauss积分公式等。2.1 梯形公式和Simpson公式梯形公式和Simpson公式是插值型积分公式。梯形公式是用一次函数逼近被积函数的积分公式，而Simpson公式是利用二次函数逼近被积函数

29、的积分公式。梯形积分公式：。计算TRAPE.F90module INTEGRAL implicit none real, parameter : PI=3.14159contains! 产生数列 subroutine GenerateData(datas, width, func) real datas(:), widthreal, external : funcreal rinteger i,nn = size(datas,1)width = PI/(n-1)r = 0do i=1,n datas(i) = func(r) r = r+widthend do end subroutine!

30、梯形法积分 real function Trape_Integral(datas, width) implicit none real datas(:) real width ! 每条数据的间隔 real SUM ! 计算所有上底加下底除以二的和 integer i,n n = size(datas,1) SUM = (datas(1)+datas(n)/2.0 do i=2,n-1 SUM=SUM+datas(i) ! 累加边长 end do Trape_Integral=SUM*width ! 计算面积和 return end function end module! 梯形法积分范例pro

31、gram main use INTEGRAL implicit none integer, parameter : N = 10 real DATAS(N), width real ANS ! 答案 real, intrinsic : sin ! 模拟用来产生数据的函数 call GenerateData(DATAS, width, sin) ANS = Trape_Integral(DATAS, width) ! 计算积分 write(*,"('ans=',F5.2)") ANS ! 写出答案 stopend programSimpson 积分公式：Sim

32、pson.F90module INTEGRAL implicit none real, parameter : PI=3.14159contains! 产生数列 subroutine GenerateData(datas, width, func) real datas(:), widthreal, external : funcreal rinteger i,nn = size(datas,1)width = PI/(n-1)r = 0do i=1,n datas(i) = func(r) r = r+widthend do end subroutine! 梯形法积分real functio

33、n Simpson_Integral(datas, width) IMPLICIT NONE real datas(:), width real sum integer i,n n = size(datas,1) if ( mod(n,2)=0 ) then write(*,*) "要有奇数条数据"stop end if sum = datas(1) + datas(n) ! 先算出头尾的和 do i=2,n-1 if ( mod(i,2)=0 ) then sum = sum + 4*datas(i) ! 把4*f(x)的部分累加起来else sum = sum + 2*

34、datas(i) ! 把2*f(x)的部分累加起来end if end do Simpson_Integral = sum * width/3.0 ! SUM再乘上H/3 就好了 returnend functionend module! SIMPSON法积分范例program main use INTEGRAL implicit none integer, parameter : N = 9 real, intrinsic : sin real datas(N), width call GenerateData(datas, width, sin) write(*,"('a

35、ns=',F6.2)") simpson_integral(datas, width) stopend program3 插值与拟合3.1 Lagrange Interpolation所谓插值，就是给定一组点和函数值，求一个函数，使得。插值一般表示成基函数的形式：，其中成作插值基函数，其满足，。Lagrange插值的基函数为：给定正弦函数的10个值，采用拉格朗日插值法求100个插值节点的值。LAGRANGE.F90module INTERPOLATE_UTILITY implicit none type point real x,y end type real, parame

36、ter : PI=3.14159 real, parameter : xmin = 0.0, xmax = PI*3.0 integer, parameter : N = 10, NP = 100 type(point) : datas(N) type(point) : interpolate(NP)contains! 产生数列 subroutine GenerateData(func)real, external : funcreal r, widthinteger iwidth = (xmax-xmin)/(N-1)r = 0do i=1,N datas(i)%x = r datas(i)

37、%y = func(r) r = r+widthend do end subroutine real function lagrange(x)real x real coeffinteger i,jlagrange = 0 do i=1,n coeff = 1 do j=1,n if ( i/=j ) coeff = coeff * (x-datas(j)%x)/(datas(i)%x-datas(j)%x) end do lagrange = lagrange + coeff*datas(i)%yend do end functionend moduleprogram main use IN

38、TERPOLATE_UTILITY implicit none real, intrinsic : sin real xinc,x integer i call GenerateData(sin) ! 产生数据点 x=0 xinc = (xmax-xmin)/(NP-1) do i=1,NP interpolate(i)%x = x interpolate(i)%y = lagrange(x) ! 插值出f(x)的值x = x+xinc end dowrite(*,"(4X'X',5X)") do I=1,Np write(*,"(F5.3,F10

39、.5)") X(I),Y(I) end do stopend program3.2 最小二乘法（Least-Square Method）最小二乘法用于数据处理，利用给定的数据拟合出函数曲线。与插值法不同的是，拟合曲线一般不通过给定的数据点，数据点分布在曲线的附近。曲线拟合根据事先假定的函数关系不同可以分为：线性和非线性拟合。按照变量的多少可以分为一元和多元拟合。假设一个铜条在0度时长1米，由于温度的不同热胀冷缩的长度变化如下：温度T（C）长度Y（米）51.047101.112151.152201.191251.252假设铜的长度关于温度是线性膨胀的，求膨胀系数。根据已知条件，可以得到

40、长度和温度的关系式为。利用最小二乘法的思想是：使得数据点到直线的距离最小。为了便于数学处理，采用距离平方和最小。N个点到直线的距离平方和为：要使距离平方和最小，将上式分别对求偏导数，并令其等于零，求出。上式等价于写成矩阵的形式为求解方程组，得到系数。程序实现如下：LEAST_SQUARE.F90module datas implicit none integer, parameter : N=5 real : temperature(N) = (/5.0,10.0,15.0,20.0,25.0/) ! 记录温度 real : length(N) = (/1.047,1.112,1.1152,1

41、.191,1.252/)! 记录不同温度下的长度 real, save : A,B ! 函数L=At+B的系数end module! 最小方差法范例program main use datas implicit none ! 调用最小方差法的子程序来计算系数A,B call least_square(temperature, length, N, A, B) write(*,"('A=',F6.4,' B=',F6.4)") A,B stopENDsubroutine least_square(x, y, n, s, t) implicit

42、none integer n real x(n) ! x上的数据点 real y(n) ! y上的数据点 real s,t ! 所要计算的系数 real A,B,C,D,E,F ! 联立方程式中的系数 integer I ! 循环计数器 ! 解 As+Bt=E 中的未知数s,t ! Cs+Dt=F ! 先设定好A,B,C,D,E,F的系数值 A = 0; B = 0; E = 0; F = 0 do I=1,N A=A+X(I)*X(I) B=B+X(I) E=E+X(I)*Y(I) F=F+Y(I) END DO C=B D=N ! 二元一次方程式有公式可解 S=(B*F-E*D)/(B*C

43、-A*D) T=(E*C-A*F)/(B*C-A*D) returnend subroutine程序执行结果：4 数据结构与算法4.1 排序4.1.1 冒泡排序法冒泡排序法是最简单的排序方法，他的步骤如下：（1）从第一个数字开始，依次把两个相邻的数据互相比较大小。如果前一个数字比后一个数字大，就交换他们的位置；（2）一直做到每一对相邻数字都比较过后，结束这一轮工作；（3）回到第一步，再作下一循环的比较。如果有n个数字要排序，就需要重复n-1次的扫描工作。冒泡排序法可以想象成让重的东西向下沉，轻的东西向上浮。等到状态稳定，就得到排序的结果。下面是冒泡排序法的程序实现：! 冒泡排序法范例! By

44、Perng 1997/8/29program BUBBLE_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 BUBBLE_SORT(A,N) ! 调用排序的子程序 write(*,"('Sort=>',10I3)") A stopend programsubroutine BUBBLE_SORT

45、(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个数据中最小的一个，把它和数列的第一个数字

46、交换位置；（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 c

47、all 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(

48、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 :

49、 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(*,"(

50、'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 ! 找到了, 返回数字在类型中的位置 SE

51、QUENTIAL_SEARCH=I returnend if end do ! 没找到时返回-1 SEQUENTIAL_SEARCH=0 returnend function4.2.2 二元搜索二元搜索法必须配合排序好的数据才能使用，假设所要寻找的数据值为K，数据存放在数组中，搜索的步骤为：（1）取出数组的中间值M和K来比较，如果M=K，结束。如果K>M，因为数组早已作好排序，所以K值一定是在数组的后半段。如果K<M，那么K值一定在数组的前半段；（2）根据K值在数组的前半段或后半段来重新分组，再回到第一步进行搜索。分组一直细分到不能再细分为止，而此时若还没有找到K，代表K值不存在。二元搜索法充分利用数组排序的特性，相比顺序搜索法可以减少搜索次数。二元搜索法的程序如下：! 折半查找法范例! By Perng 1997/8/31program BINARY_SE