常微分方程初值問題數(shù)值解法

1、常微分方程初值問題數(shù)值解法朱欲輝（浙江海洋學院 數(shù)理信息學院, 浙江 舟山 316004）摘要：在常微分方程的課程中討論的都是對一些典型方程求解析解的方法. 然而在生產(chǎn)實際和科學研究中所遇到的問題往往很復雜, 在很多情況下都不可能給出解的解析表達式. 本篇文章詳細介紹了常微分方程初值問題的一些數(shù)值方法, 導出了若干種數(shù)值方法, 如Euler法、改進的Euler法、RungeKutta法以及線性多步法中的Adams顯隱式公式和預測校正公式, 并且對其穩(wěn)定性及收斂性作了理論分析. 最后給出了數(shù)值例子, 分別用不同的方法計算出近似解, 從得出的結果對比各種方法的優(yōu)缺點. 關鍵詞：常微分方程; 初值問

2、題; 數(shù)值方法; 收斂性; 穩(wěn)定性; 誤差估計Numerical Method for Initial-Value ProblemsZhu Yuhui(School of Mathematics, Physics, and Information Science, Zhejiang Ocean University, Zhoushan, Zhejiang 316004)Abstract: In the course about ordinary differential equations, the methods for analytic solutions of some typical

3、equations are often discussed. However, in scientific research, the problems are very complex and the analytic solutions about these problems cant be expressed explicitly. In this paper, some numerical methods for the initial-value problems are introduced. these methods include Euler method, improve

4、d Euler method, Runge-Kutta method and some linear multistep method (e.g. Adams formula and predicted-corrected formula). The stability and convergence about the methods are presented. Some numerical examples are give to demonstrate the effectiveness and accuracy of theoretical analysis. Keywords: O

5、rdinary differential equation; Initial-value problem; Numerical method; Convergence; Stability; Error estimate1前言自然界和工程技術中的很多現(xiàn)象, 例如自動控制系統(tǒng)的運行、電力系統(tǒng)的運行、飛行器的運動、化學反應的過程、生態(tài)平衡的某些問題等, 都可以抽象成為一個常微分方程初值問題. 其真解通常難以通過解析的方法來獲得, 至今有許多類型的微分方程還不能給出解的解析表達式, 一般只能用數(shù)值的方法進行計算. 有關這一問題的研究早在十八世紀就已經(jīng)開始了, 特別是計算機的普遍應用, 許多微分方程問

6、題都獲得了數(shù)值解, 從而能使人們認識解的種種性質及其數(shù)值特征, 為工程技術等實際問題提供了定量的依據(jù).關于常微分方程初值問題的數(shù)值計算方法, 許多學者己經(jīng)做了大量的工作. 1768年, Euler提出了關于常微分方程初值問題的方法, 1840年, Cauchy第一次對初值問題進行了仔細的分析, 早期的常微分方程數(shù)值解的問題來源于天體力學. 在1846年, 當Adams還是一個學生的時候, 和Le Verrier一起根據(jù)天王星軌道中出現(xiàn)的己知位置, 預測了它下一次出現(xiàn)的位置. 1883年, Adams提出了Adams一Bashforth和Adams一Moulton方法. Rull(1895年)、

7、Heun(1900年)和Kutta(1901年)提出Runge.Kutta方法.二十世紀五十年代, Dahlquist建立了常微分方程數(shù)值解法的穩(wěn)定性理論, 線性多步法是常微分方程初值問題的一種數(shù)值方法. 由于通常的數(shù)值方法, 其絕對穩(wěn)定區(qū)域是有限的, 不適用于求解剛性常微分的初值問題. 剛性微分方程常常出現(xiàn)于航空、航天、熱核反應、自動控制、電子網(wǎng)絡及化學動力學等一系列與國防和現(xiàn)代化建設密切相關的高科技領域, 具有無容置疑的重要性. 因此, 剛性微分方程的研究工作早在二十世紀五十年代就開始了, 1965年, 在愛丁堡舉行的IFIP會議后, 更進一步地認識剛性方程的普遍性和重要性. 自從六十年代

8、初, 許多數(shù)值分析家致力于探討剛性問題的數(shù)值方法及其理論, 注意到剛性問題對傳統(tǒng)數(shù)值積分方法所帶來的挑戰(zhàn). 這一時期, 人們的研究主要集中在算法的線性穩(wěn)定性上, 就是基于試驗方程數(shù)值解的穩(wěn)定性研究. 在此領域發(fā)表了大量的論文, 取得了許多重要的理論成果. 例如, 1963年, Dahlquist給出A穩(wěn)定性理論, 1967年, Widlund給出穩(wěn)定性理論, 1969年, Gear將穩(wěn)定性減弱, 給出剛性(Stiff)穩(wěn)定性理論, 并找到了當?shù)膋步k階的剛性穩(wěn)定方法, 1969年Dill找到剛性穩(wěn)定的7階和8階以及1970年Jain找到剛性穩(wěn)定的9階到11階, 但可用性沒有檢驗. 這些穩(wěn)定性理

9、論和概念都是在線性試驗方程的框架下推導出的, 從嚴格的數(shù)學意義上來說, 這些理論只適用于常系數(shù)線性自治系統(tǒng). 但從實用的觀點來說, 這些理論無疑是合理和必要的, 對剛性問題的算法設計具有重要的指導意義. 在八十至九十年代, 國內也有一些學者研究線性理論, 主要有匡蛟勛、陳果良、項家祥、李壽佛、黃乘明、李慶揚和費景高等.線性理論雖然對一般問題具有指導作用, 但其不能作為非線性剛性問題算法的穩(wěn)定性理論研究基礎. 為了將線性理論推廣到非線性問題中, 人們開始對非線性模型問題進行研究.但是, 早期文獻主要致力于數(shù)值方法基于經(jīng)典Lipschitz條件下的經(jīng)典收斂理論, 即認為良好的穩(wěn)定性加上經(jīng)典相容性和

10、經(jīng)典相容階就足以描述方法的整體誤差性態(tài). 直到1974 年, Prothero和Robinson首先注意到算法的經(jīng)典誤差估計由于受剛性問題巨大參數(shù)的影響而嚴重失真, 產(chǎn)生階降低現(xiàn)象, 這時人們認識到經(jīng)典收斂理論對于非線性剛性問題以及線性模型的不足. 于是, 1975年, Dahlquist和Butcher分別提出了單支方法和線性多步法的G一穩(wěn)定概念和B穩(wěn)定概念. 這兩個概念填補了非線性穩(wěn)定性分析理論, 引起了計算數(shù)學家們的極大關注, 在上述理論的基礎上, 1975年至1979年, Burrage和Butcher提出了AN一穩(wěn)定性與BN一穩(wěn)定性概念, 并相應地建立了基本的B一穩(wěn)定及代數(shù)穩(wěn)定理論.

11、 1981至1985年, Frank, schneid和ueberhuber建立Runge一kutta方法的B一收斂理論. B穩(wěn)定與B一收斂理論統(tǒng)稱B一理論, 它是常微分方程數(shù)值解法研究領域的巨大成就之一, 是剛性問題算法理論的突破性進展, 標志著剛性問題研究從線性向非線性情形深入發(fā)展. 國內也有眾多學者致力于B一理論的研究, 如李壽佛、曹學年等. 1989年, 李壽佛將Dahlquist的G一穩(wěn)定概念推廣到更一般的(C,P,Q)代數(shù)穩(wěn)定, 克服了G穩(wěn)定的線性多步法不能超過二階的限制. 對于一般線性方法, 李壽佛建立了一般線性方法的(K,P,Q)穩(wěn)定性理論及(K,P,Q)弱代數(shù)穩(wěn)定準則和多步R

12、unge-Kutta法的一系列代數(shù)準則. 此外, Dahquist, Butcher和Hairer分別深刻地揭示了單支方法、一般線性方法和Runge.Kutta方法線性與非線性穩(wěn)定性之間的內在關系. 為了求 解剛性微分方程, 不少文獻中構造含有穩(wěn)定參數(shù)的線性多步方法, 利用適當選擇穩(wěn)定參數(shù)來擴大方法的穩(wěn)定區(qū)域. 所有改進的思想, 都是通過構造一些特殊的顯式或隱式線性多步法, 使其具有增大的穩(wěn)定域, 或使一穩(wěn)定的角增大. 八十年代, 就成為國內外學者所研究的一個課題, 學者主要有Rodabaugh和Thompson、Feinberg、李旺堯、李壽佛、包雪松、徐洪義、劉發(fā)旺、匡蛟勛、項家祥、蔣立新

13、、李慶揚、謝敬東和李林忠等.當前國內外研究剛性問題的一個主要趨勢就是在B一理論指導下尋找更為有效的新算法. 另一個發(fā)展趨勢就是力圖突破單邊lipschitz常數(shù)和內積范數(shù)的局限, 建立比B一理論更為普遍的定量分析收斂理論. 近年來, 剛性延遲系統(tǒng)的算法研究成為剛性問題的另一個熱點研究領域, 張誠堅將Burrage等人創(chuàng)立的針對剛性常微分系統(tǒng)的B一理論拓展到非線性剛性延遲系統(tǒng).常微分方程的數(shù)值算法發(fā)展到今天己有了線性多步法、龍格一庫塔法和在此基礎上發(fā)展起來的單支方法、分塊方法、循環(huán)方法、外推法、混合方法、二階導數(shù)法以及各種常用的估校正算法. 其中經(jīng)常用到的線性多步法公式有Euler公式、Heun

14、公式、中點公式、Milne公式、Adams公式、simpson公式、Hamming公式, Gear方法、Adams預估一校正法和Mile預估一Hamming校正法公式等, 此外還包含許多迄今尚末探明的新公式. Burage曾將線性多步法和RungeKutta法比作大海中的兩座小島, 在浩瀚的汪洋之中, 還有許多到現(xiàn)在沒有發(fā)現(xiàn)的新方法.本篇文章詳細介紹了常微分方程初值問題的一些數(shù)值方法, 導出了若干種數(shù)值方法, 如Euler法、改進的歐拉法、顯式龍格-庫塔法、隱式龍格-庫塔法以及線性多步法中的Adams顯隱式公式和預測校正公式, 并且對其穩(wěn)定性及收斂性作了理論分析. 最后給出了數(shù)值例子, 進行了

15、計算機程序算法的分析與實現(xiàn), 以計算機的速度優(yōu)勢來彌補計算量大的不足. 從得出的結果對比各種方法的優(yōu)缺點. 2常微分方程初值問題的數(shù)值解法2.1 數(shù)值方法的基本思想考慮一階常微分方程初值問題 (2.1)的數(shù)值解法, 其中是和的己知函數(shù), 是給定的初始值.對于常微分方程初值問題(2.1)數(shù)值解法, 就是要算出精確解在區(qū)間上的一系列離散節(jié)點的函數(shù)值, , 的近似值, , , .相鄰兩個節(jié)點的間距稱為步長, 這時節(jié)點也可以表示為. 數(shù)值解法需要把連續(xù)性的問題加以離散化, 從而求出離散節(jié)點的數(shù)值解.通常微分方程初值問題(2.1)的數(shù)值方法可以分為兩類:(1) 單步法-計算在處的值僅取決于處的應變量及其

16、導數(shù)值.(2) 多步法-計算在處的值需要應變量及其導數(shù)在之前的多個網(wǎng)個節(jié)點出的值.2.2 Euler方法2.2.1 Euler公式若將函數(shù)在點處的導數(shù)用兩點式代替, 即,再用近似地代替, 則初值問題(2.1)變?yōu)?(2.2)(2.2)式就是著名的歐拉(Euler)公式.2.2.2 梯形公式歐拉法形式簡單精度低, 為了提高精度, 對方程的兩端在區(qū)間上積分得 (2.3)改用梯形方法計算其積分項, 即代入式(2.3), 并用近似代替式中即可得到梯形公式 (2.4)由于數(shù)值積分的梯形公式比矩形公式精度高, 所以梯形公式(2.4)比歐拉公式(2.2)的精度高一個數(shù)值方法.式(2.4)的右端含有未知的,

17、它是關于的函數(shù)方程, 這類方法稱為隱式方法.2.2.3 改進的歐拉公式梯形公式實際計算時要進行多次迭代, 因而計算量較大. 在實用上, 對于梯形公式(2.4)只迭代一次, 即先用歐拉公式算出的預估值, 再用梯形公式(2.4)進行一次迭代得到校正值, 即 (2.5)2.2.4 歐拉法的局部截斷誤差衡量求解公式好壞的一個主要標準是求解公式的精度, 因此引入局部截斷誤差和階數(shù)概念.定義2.1 在準確的前提下, 即時, 用數(shù)值方法計算的誤差, 稱為該數(shù)值方法計算時的局部截斷誤差.對于歐拉公式, 假定, 則有而將真解在處按二階泰勒展開式有因此有定義2.2 若數(shù)值方法的局部截斷誤差為, 則稱這種數(shù)值方法的

18、階數(shù)是P. 步長(h<1)越小, P越高, 則局部截斷誤差越小, 計算精度越高.2.3 RungeKutta方法2.3.1 RungeKutta方法的基本思想歐拉公式可改寫成則的表達式與的泰勒展開式的前兩項完全相同, 即局部截斷誤差為.改進的歐拉公式又可改寫成.上述兩組公式在形式上有一個共同點: 都是用在某些點上值的線形組合得出的近似值, 而且增加計算的次數(shù), 可提高截斷誤差的階. 如歐拉公式, 每步計算一次的值, 為一階方法. 改進的歐拉公式需計算兩次的值, 它是二階方法. 它的局部截斷誤差為.于是可考慮用函數(shù)在若干點上的函數(shù)值的線形組合來構造近似公式, 構造時要求近似公式在處的泰勒展

19、開式與解在處的泰勒展開式的前面幾項重合, 從而使近似公式達到所需要的階數(shù). 既避免求偏導, 又提高了計算方法精度的階數(shù). 或者說, 在這一步內多預報幾個點的斜率值, 然后將其加權平均作為平均斜率. 則可構造出更高精度的計算格式, 這就是RungeKutta方法的基本思想.2.3.2 RungeKutta方法的構造一般地, RungeKutta方法設近似公式為 （下面的公式修改了） 其中, , 都是參數(shù), 確定它們的原則是使近似公式在處的泰勒展開式與在處的泰勒展開式的前面的項盡可能多的重合, 這樣就使近似公式有盡可能高的精度. 以此我們可以通過一個復雜的計算過程得出常用的的三階和四階RungeK

20、utta公式和 (2.6)式(2.6)稱為經(jīng)典RungeKutta方法.2.4 線性多步法在逐步推進的求解過程中, 計算之前事實上已經(jīng)求出了一系列的近似值, , , . 如果充分利用前面多步的信息來預測, 則可期望獲得較高的精度, 這就是構造多步法的基本思想.線性k步方法的一般公式為 (2.7)其中, 均為與n無關的常數(shù), . 當時為顯格式; 當時為隱格式. 特別當時為Euler公式；當時為梯形公式.定義 2.3 稱為k步公式(2.7)在處的局部截斷誤差. 當時稱式(2.7)是p階的.應用方程可知局部截斷誤差也可寫成為定義2.4 如果線性k步方法(2.7)至少是1階的, 則稱是相容的; 如果線

21、性k步法(2.7)是p階的, 則稱是p階相容的.2.4.1 Adams外插法將微分方程的兩端從到進行積分, 得到 (2.8)我們用插值多項式代替右端的被積函數(shù).Adams外插法選取k個點作為插值基點構造的k-1階多項式Adams外插法的計算公式為其中滿足如下代數(shù)遞推式:,根據(jù)此遞推公式, 可逐個的計算, 表2.1給出了的部分數(shù)值:表2.1j01234511/25/123/8251/72095/2882.4.2 Adams內插法根據(jù)插值理論知道, 插值節(jié)點的選擇直接影響著插值公式的精度, 同樣次數(shù)的內插公式的精度要比外插公式的高.仍假定已按某種方法求得問題(2.1)的解在處的數(shù)值, 并選取插值節(jié)

22、點, p是正整數(shù), 用Lagrange型插值多項式構造可以導出解初值問題(2.1)的Adams內插公式: (2.9)當時上式就退化為內插公式.內插公式(2.9)除了包含在處的已知值外, 還包含了在點, 處的未知值.因此內插公式(2.9)只給出了未知量的關系式, 實際計算時仍需要解聯(lián)立方程組. 的內插公式是最適用的, 采用Newton向后插值公式得到Adams內插公式其中系數(shù)定義為其中滿足如下代數(shù)遞推式:根據(jù)此遞推公式, 可逐個的計算, 表2.2給出了的部分數(shù)值:表2.2j0123451-1/2-1/12-1/24-19/720-3/1602.5 算法的收斂性和穩(wěn)定性2.5.1 相容性初值問題(

23、2.1):的顯式單步法的一般形式為, (2.10)其中, . 這樣我們用差分方程初值問題(2.10)的解作為問題(2.1)的解在處的近似值, 即.因此, 只有在問題(2.1)的解使得逼近時, 才有可能使(2.10)的解逼近于問題（2.1）的解. 從而, 我們期望對任一固定的 都有假設增量函數(shù)關于h連續(xù), 則有定義 2.5 若關系式成立, 則稱單步法(2.10)與微分方程初值問題(2.1)相容.容易驗證, Euler法與改進Euler法均滿足相容性條件. 事實上, 對Euler法, 增量函數(shù)為自然滿足相容性條件, 改進Euler法的增量函數(shù)為因為, 從而有所以Euler法與改進Euler法均與初

24、值問題（2.1）相容. 一般的, 如果顯示單步法有p階精度, 則其局部誤差為上式兩端同除以h, 得令, 如果, 則有即所以的顯示單步法均與初值問題(2.1)相容.由此各階RK方法與初值問題(2.1)相容.2.5.2 收斂性定義2.6 一種數(shù)值方法稱為是收斂的, 如果對于任意初值及任意 都有 其中為初值問題(2.1)的準確解.按定義, 數(shù)值方法的收斂性需根據(jù)該方法的整體截斷誤差來判定. 已知Euler方法的整體截斷誤差有估計式當, , 故Euler方法收斂.定理 2.1 設顯式單步法具有p階精度, 其增量函數(shù)關于y滿足Lipschitz條件, 問題(2.1)是精確的, 既, 則顯式單步法的整體截

25、斷誤差為.定理 2.2設增量函數(shù)在區(qū)域S上連續(xù), 且關于y滿足Lipschitz條件時, 則顯式單步法收斂的充分必要條件是相容性條件成立, 即.2.5.3 穩(wěn)定性定義 2.7 如果存在正常數(shù)及C, 使得對任意的初始出發(fā)值, 單步法(2.10)的相應精確解, 對所有的, 恒有, 則稱單步法是穩(wěn)定的.定理2.3若對于, 以及一切實數(shù)y, 關于y滿足Lipschitz條件, 則單步法(2.7)是穩(wěn)定的.定義 2.8對給定的微分方程和給定的步長h, 如果由單步法計算時有大小為的誤差, 即計算得出, 而引起其后之值的變化小于, 則說該單步法是絕對穩(wěn)定的.一般只限于典型微分方程考慮數(shù)值方法的絕對穩(wěn)定性,

26、其中為復常數(shù)(我們僅限于為實數(shù)的情形). 若對于所有的, 單步法都絕對穩(wěn)定, 則稱為絕對穩(wěn)定區(qū)間.根據(jù)以上定義我們可以得出Euler方法的絕對穩(wěn)定區(qū)間為(-2,0), 梯形公式的絕對穩(wěn)定區(qū)間為, RungeKutta的絕對穩(wěn)定區(qū)間為(-2.78,0).3 實例分析例3.1分別用Euler法、改進Euler法、RungeKutta解初值問題的數(shù)值解, 取h=0.1, N=10, 并與真實解作出比較.由于該簡單方程可以用數(shù)學方法求得其精確描述式, 所以可以據(jù)此檢驗近似數(shù)值解同真實解的誤差情況. 對于其他一些結構復雜的常微分方程的數(shù)值解實現(xiàn)方法也是一樣的.（1）使用歐拉算法進行一般求解經(jīng)過歐拉算法程

27、序計算得出在各點處的近似數(shù)值解及各自同精確解的誤差, 如下表3.1表3.1各分點(數(shù)值解)(精確解)(誤差)0.11.0000000.9900990.0099010.20.9800000.9615380.0184620.30.9415840.9174310.0241530.40.8883890.8620690.0263200.50.8252500.8000000.0252500.60.7571470.7352940.0218520.70.6883540.6711410.0172130.80.6220180.6097560.0122620.90.5601130.5524860.0076261.0

28、0.5036420.5000000.0036421.10.4529110.4524890.0004221.20.4077830.4098360.002053從實驗結果看誤差已經(jīng)不算太小, 況且這還僅僅是一個用于實驗的比較簡單的常微分方程, 對于實際工程中應用的結構復雜的方程其求解結果的誤差要遠比此大得多, 由于還存在著局部截斷誤差和整體截斷誤差, 因此可采取一定措施來抑制減少誤差, 盡量使結果精確.（2）使用改進歐拉算法進行一般求解經(jīng)過改進歐拉算法程序計算得出在各點處的近似數(shù)值解及各自同精確解的誤差, 如下表3.2表3.2各分點(數(shù)值解)(精確解)(誤差)0.10.9900000.990099

29、0.0000990.20.9613660.9615380.0001730.30.9172460.9174310.0001850.40.8619540.8620690.0001150.50.8000340.8000000.0000340.60.7355270.7352940.0002330.70.6715870.6711410.0004460.80.6103990.6097560.0006430.90.5532890.5524860.0008031.00.5009190.5000000.0009191.10.4534790.4524890.0009901.20.4108590.4098360.0

30、01023可以看出, 這種經(jīng)過改進的歐拉算法所存在的誤差已大為減少, 誤差的減少主要是由于先利用了歐拉公式對的值進行了預估, 然后又利用梯形公式對預估值作了校正, 從而在預估校正的過程中減少了誤差. 此算法有一定的優(yōu)點, 在一些實際而且簡單的工程計算中可以直接單獨應用.（3）使用經(jīng)典RungeKutta法進行一般求解經(jīng)過經(jīng)典RungeKutta法程序計算得出在各點處的近似數(shù)值解及各自同精確解的誤差, 如下表3.3表3.3各分點(數(shù)值解)(精確解)(誤差)0.10.9900990.9900990.0000000.20.9615380.9615380.0000000.30.9174310.9174

31、310.0000010.40.8620680.8620690.0000010.50.7999990.8000000.0000010.60.7352940.7352940.0000010.70.6711410.6711410.0000000.80.6097560.6097560.0000000.90.5524870.5524860.0000001.00.5000010.5000000.0000011.10.4524890.4524890.0000011.20.4098370.4098360.000001從表3.3記錄的程序運行結果來看, 所計算出來的常微分方程的數(shù)值解是非常精確的, 用浮點數(shù)進行

32、運算時由近似所引入的誤差幾乎不會對計算結果產(chǎn)生影響, 這樣高的精度是很令人滿意的. 無論是從計算速度上還是從計算精度要求上, 都能很好地滿足工程計算的需要.例3.2分別用Euler法、改進Euler法、RungeKutta解初值問題的數(shù)值解, 取等分數(shù)分別為N=12, 24, 36, ,120, 分別計算出與真實解的最大誤差.（1） 使用歐拉算法進行一般求解經(jīng)過計算我們可以得出初值問題在取不通步長時, 數(shù)值解與真實解的最大誤差如表3.4表3.4N(等分數(shù))最大誤差N(等分數(shù))最大誤差120.026320720.004052240.012547840.003465360.008233960.00

33、3027480.0061261080.002687600.0048781200.002416從表3.4記錄的程序運行結果來看當?shù)确謹?shù)N變大時, 它的誤差正在減小, 根據(jù)定義2.6我們可以證得該方法是收斂的.（2）使用改進歐拉算法進行一般求解經(jīng)過計算我們可以得出初值問題在取不通步長時, 數(shù)值解與真實解的最大誤差如表3.5表3.5N(等分數(shù))最大誤差N(等分數(shù))最大誤差120.001023720.000028240.000255840.000021360.000113960.000016480.0000631080.000013600.0000411200.000010從表3.5記錄的程序運行結果

34、來看當?shù)确謹?shù)N變大時, 它的誤差正在減小, 根據(jù)定義2.6我們可以證得該方法是收斂的.（3）使用經(jīng)典RungeKutta法進行一般求解經(jīng)過計算我們可以得出初值問題在取不通步長時, 數(shù)值解與真實解的最大誤差如表3.6表3.6N(等分數(shù))最大誤差N(等分數(shù))最大誤差120.000001720.000000240.000000840.000000360.000000960.000000480.0000001080.000000600.0000001200.000000從表3.6記錄的程序運行結果來看當?shù)确謹?shù)N變大時, 它的誤差正在減小, 根據(jù)定義2.6我們可以證得該方法是收斂的.（4）收斂速度對比為對比以上三種方法的收斂速度, 我們計算出了它們的誤差精度達到的最小等分數(shù)N如下表3.7表3.7方法N(最小等分數(shù))誤差精度Euler法2900改進Euler法40RungeKutta法5從表3.7記錄的程序運行結果來看RungeKutta法的收斂速度最快, 改進Euler法其次, 而Euler法最差. 由此看來RungeKutta法是他們當中最理想的數(shù)值解法.小結針對工程計算中遇到的常微分方程初值問題的求解, 根據(jù)實際情況引如了計算機作為輔助計算工具, 并根據(jù)高等數(shù)學的有關知識將歐拉公式、改進的歐拉公式、經(jīng)典龍格庫塔方法等融入到計算機程序算法