Fortran基础编程——牛顿拉弗森迭代

微信公众号：[易木木响叮当]
关注可了解更多的[有限元数值仿真技巧]。问题或建议，请公众号留言;
如果你觉得木木同学对你有帮助，欢迎赞赏。

今日推文主要分享一个非线性方程的牛顿-拉弗森迭代解法，借助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也是同样的道理。

1. use Newton:主程序中需要声明的变量存储在Module NEWTON中，使用的时候只需use即可；

2. OPEN(UNIT=11,FILE='FOUT1510.TXT')：表示在当前目录里面打开FOUT1510.TXT文件，若没有则表示新建一个FOUT1510.TXT文件，记通道号为11。
   关于通道号的解释，书中有相关讲解，木木在学习的过程中，将之理解为与write语句有关，打开的通道号用于写入或者读入该通道，比如WRITE(11,46)X,ITER表示按照48行代码标号的格式写入到11号通道文件中，可以自己改动一下代码try one try；

3. CALL SOLVE(X,ITER)表示调用牛顿-拉弗森迭代子程序；

4. 46    FORMAT(T5,'Newton迭代法计算方程的根',//,&：46表示行代码，多用于write语句中，表示格式输出位置，T5表示在该行的第5个字符位置以后进行输出，/表示换行，//表示换2行，行尾&表示续行；

5. 3X表示空出3个字符位置，F15.10表示实数字符长度为15个字符，小数点后面的占10位，I5表示整数位置为5个字符。

Moudle

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”值

今日有关Fortran编程的知识就在这里，完整的源文件可在后台回复：Fortran_Newton我们下次再见，欢迎转发与小赞~

最后插一句话，木木最近参加了技术邻的仿真比赛，需要小伙伴的投票支持，大家有时间的话可以动动小手点点投票（电脑端打开更加顺畅，手机端可能有Bug）。等木木获奖后，再为大家继续抽奖活动，更多福利等着大家~

投票链接：基于扩展有限元的混凝土受力开裂分析