今日推文主要分享一个非线性方程的牛顿-拉弗森迭代解法,借助Fortran语言,讲述Fortran编程时需要注意的地方。理论及在Abaqus中的实现过程已在上几期推文基于Abaqus的Newton-Raphson算法中说明,本次主要说明Fortran编程时需要注意的地方,本文代码主要参考:《Fortran程序设计权威指南》。
例:
主程序:Solve.f90
,子程序:New_Raphson.f90
,函数文件:function.f90
,Module模块:NEWTON.F90
。将整个程序分块编写,在主程序中调用即可,大型的Fortran中显得尤为重要。
主程序代码:
PROGRAM main
! 主程序:Newton迭代法计算方程的根
use NEWTON
OPEN(UNIT=11,FILE='FOUT1510.TXT')
OPEN(UNIT=12,FILE='IM_RESULT1510.TXT')
CALL SOLVE(X,ITER)
WRITE(11,46)X,ITER
46 FORMAT(T5,'Newton迭代法计算方程的根',//,&
3X,'X= ',F15.10,/,&
3X,'ITER=',I5)
END PROGRAM main
由上述程序可看出,主程序只有use、call
、输出语句,甚至可以将输出语句编一个子程序,在主程序中用call
调用即可,如此以来,使得自己的代码更加简洁明朗,在编写Python、Matlab也是同样的道理。
use Newton
:主程序中需要声明的变量存储在Module NEWTON
中,使用的时候只需use
即可;
OPEN(UNIT=11,FILE='FOUT1510.TXT')
:表示在当前目录里面打开FOUT1510.TXT文件,若没有则表示新建一个FOUT1510.TXT文件,记通道号为11。
关于通道号的解释,书中有相关讲解,木木在学习的过程中,将之理解为与write
语句有关,打开的通道号用于写入或者读入该通道,比如WRITE(11,46)X,ITER
表示按照48行代码标号的格式写入到11号通道文件中,可以自己改动一下代码try one try;
CALL SOLVE(X,ITER)
表示调用牛顿-拉弗森迭代子程序;
46 FORMAT(T5,'Newton迭代法计算方程的根',//,&
:46表示行代码,多用于write
语句中,表示格式输出位置,T5
表示在该行的第5个字符位置以后进行输出,/
表示换行,//
表示换2行,行尾&
表示续行;
3X
表示空出3个字符位置,F15.10
表示实数字符长度为15个字符,小数点后面的占10位,I5
表示整数位置为5个字符。
MODULE NEWTON
IMPLICIT REAL*8(A-Z)
INTEGER ITER
END MODULE NEWTON
Module中可放置程序所需要声明的所有变量,也可以是子程序、函数,用的时候仅需use
语句调用即可,不只主程序可以使用,子程序中也可以使用,是fortran90以后添加的功能,大大提升编程的效率。
SUBROUTINE SOLVE(X,ITER)
! Newton迭代法求方程的根
IMPLICIT REAL*8(A-Z)
INTEGER I,ITER
INTEGER:: MAX=200
PRINT*,"x0=",x0,"tol=",tol
X1=X0
DO I=1,MAX
X2=X1-FUNC(X1)/DFUNC(X1)
DX=DABS(X2-X1)
! 输出中间结果
WRITE(12,*)I,X2,FUNC(X2)
IF (DX<TOL) EXIT
X1=X2
END DO
X=X2
ITER=I
END SUBROUTINE SOLVE
具体的细节看对照程序捋一捋,木木在此讲一个有趣的现象:在源程序没有指定初值 print
出来,结果发现:x0=4.2439915819305446E-314
tol=6.3659873728958169E-314
原来系统在没有指定变量初值时,会指定一个近似于“0”值