微信公眾號(hào)：[易木木響叮當(dāng)]
關(guān)注可了解更多的[有限元數(shù)值仿真技巧]。問題或建議，請(qǐng)公眾號(hào)留言;
如果你覺得木木同學(xué)對(duì)你有幫助，歡迎贊賞。

今日推文主要分享一個(gè)非線性方程的牛頓-拉弗森迭代解法，借助Fortran語言，講述Fortran編程時(shí)需要注意的地方。理論及在Abaqus中的實(shí)現(xiàn)過程已在上幾期推文基于Abaqus的Newton-Raphson算法中說明，本次主要說明Fortran編程時(shí)需要注意的地方，本文代碼主要參考：《Fortran程序設(shè)計(jì)權(quán)威指南》。

 例：

主程序

主程序：Solve.f90，子程序：New_Raphson.f90，函數(shù)文件：function.f90，Module模塊：NEWTON.F90。將整個(gè)程序分塊編寫，在主程序中調(diào)用即可，大型的Fortran中顯得尤為重要。
主程序代碼：

PROGRAM main
    ! 主程序：Newton迭代法計(jì)算方程的根
    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迭代法計(jì)算方程的根',//,&
        3X,'X=   ',F15.10,/,&
        3X,'ITER=',I5)
END PROGRAM main

由上述程序可看出，主程序只有use、call、輸出語句，甚至可以將輸出語句編一個(gè)子程序，在主程序中用call調(diào)用即可，如此以來，使得自己的代碼更加簡潔明朗，在編寫Python、Matlab也是同樣的道理。

1. use Newton:主程序中需要聲明的變量存儲(chǔ)在Module NEWTON中，使用的時(shí)候只需use即可；

2. OPEN(UNIT=11,FILE='FOUT1510.TXT')：表示在當(dāng)前目錄里面打開FOUT1510.TXT文件，若沒有則表示新建一個(gè)FOUT1510.TXT文件，記通道號(hào)為11。
   關(guān)于通道號(hào)的解釋，書中有相關(guān)講解，木木在學(xué)習(xí)的過程中，將之理解為與write語句有關(guān)，打開的通道號(hào)用于寫入或者讀入該通道，比如WRITE(11,46)X,ITER表示按照48行代碼標(biāo)號(hào)的格式寫入到11號(hào)通道文件中，可以自己改動(dòng)一下代碼try one try；

3. CALL SOLVE(X,ITER)表示調(diào)用牛頓-拉弗森迭代子程序；

4. 46    FORMAT(T5,'Newton迭代法計(jì)算方程的根',//,&：46表示行代碼，多用于write語句中，表示格式輸出位置，T5表示在該行的第5個(gè)字符位置以后進(jìn)行輸出，/表示換行，//表示換2行，行尾&表示續(xù)行；

5. 3X表示空出3個(gè)字符位置，F15.10表示實(shí)數(shù)字符長度為15個(gè)字符，小數(shù)點(diǎn)后面的占10位，I5表示整數(shù)位置為5個(gè)字符。

Moudle

MODULE NEWTON
IMPLICIT REAL*8(A-Z)
INTEGER ITER 
END MODULE  NEWTON

Module中可放置程序所需要聲明的所有變量，也可以是子程序、函數(shù)，用的時(shí)候僅需use語句調(diào)用即可，不只主程序可以使用，子程序中也可以使用，是fortran90以后添加的功能，大大提升編程的效率。

子程序、函數(shù)

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)
        ! 輸出中間結(jié)果
        WRITE(12,*)I,X2,FUNC(X2)
        IF (DX<TOL) EXIT
        X1=X2
    END DO
    X=X2
    ITER=I
END SUBROUTINE SOLVE

具體的細(xì)節(jié)看對(duì)照程序捋一捋，木木在此講一個(gè)有趣的現(xiàn)象：在源程序沒有指定初值   和精度   ，程序卻可以跑通，木木就納悶了，系統(tǒng)怎么知道這兩個(gè)值的？于是我就將這兩個(gè)變量print出來，結(jié)果發(fā)現(xiàn)：
x0=4.2439915819305446E-314
tol=6.3659873728958169E-314
原來系統(tǒng)在沒有指定變量初值時(shí)，會(huì)指定一個(gè)近似于“0”值

今日有關(guān)Fortran編程的知識(shí)就在這里，完整的源文件可在后臺(tái)回復(fù)：Fortran_Newton我們下次再見，歡迎轉(zhuǎn)發(fā)與小贊~

最后插一句話，木木最近參加了技術(shù)鄰的仿真比賽，需要小伙伴的投票支持，大家有時(shí)間的話可以動(dòng)動(dòng)小手點(diǎn)點(diǎn)投票（電腦端打開更加順暢，手機(jī)端可能有Bug）。等木木獲獎(jiǎng)后，再為大家繼續(xù)抽獎(jiǎng)活動(dòng)，更多福利等著大家~

投票鏈接：基于擴(kuò)展有限元的混凝土受力開裂分析