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

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

2、題; 數(shù)值方法; 收斂性; 穩(wěn)定性; 誤差估計(jì)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前言自然界和工程技術(shù)中的很多現(xiàn)象, 例如自動(dòng)控制系統(tǒng)的運(yùn)行、電力系統(tǒng)的運(yùn)行、飛行器的運(yùn)動(dòng)、化學(xué)反應(yīng)的過程、生態(tài)平衡的某些問題等, 都可以抽象成為一個(gè)常微分方程初值問題. 其真解通常難以通過解析的方法來獲得, 至今有許多類型的微分方程還不能給出解的解析表達(dá)式, 一般只能用數(shù)值的方法進(jìn)行計(jì)算. 有關(guān)這一問題的研究早在十八世紀(jì)就已經(jīng)開始了, 特別是計(jì)算機(jī)的普遍應(yīng)用, 許多微分方程問

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

22、點(diǎn), p是正整數(shù), 用Lagrange型插值多項(xiàng)式構(gòu)造可以導(dǎo)出解初值問題(2.1)的Adams內(nèi)插公式: (2.9)當(dāng)時(shí)上式就退化為內(nèi)插公式.內(nèi)插公式(2.9)除了包含在處的已知值外, 還包含了在點(diǎn), 處的未知值.因此內(nèi)插公式(2.9)只給出了未知量的關(guān)系式, 實(shí)際計(jì)算時(shí)仍需要解聯(lián)立方程組. 的內(nèi)插公式是最適用的, 采用Newton向后插值公式得到Adams內(nèi)插公式其中系數(shù)定義為其中滿足如下代數(shù)遞推式:根據(jù)此遞推公式, 可逐個(gè)的計(jì)算, 表2.2給出了的部分?jǐn)?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)的解使得逼近時(shí), 才有可能使(2.10)的解逼近于問題（2.1）的解. 從而, 我們期望對(duì)任一固定的 都有假設(shè)增量函數(shù)關(guān)于h連續(xù), 則有定義 2.5 若關(guān)系式成立, 則稱單步法(2.10)與微分方程初值問題(2.1)相容.容易驗(yàn)證, Euler法與改進(jìn)Euler法均滿足相容性條件. 事實(shí)上, 對(duì)Euler法, 增量函數(shù)為自然滿足相容性條件, 改進(jìn)Euler法的增量函數(shù)為因?yàn)? 從而有所以Euler法與改進(jìn)Euler法均與初

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

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

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

27、序計(jì)算得出在各點(diǎn)處的近似數(shù)值解及各自同精確解的誤差, 如下表3.1表3.1各分點(diǎn)(數(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從實(shí)驗(yàn)結(jié)果看誤差已經(jīng)不算太小, 況且這還僅僅是一個(gè)用于實(shí)驗(yàn)的比較簡(jiǎn)單的常微分方程, 對(duì)于實(shí)際工程中應(yīng)用的結(jié)構(gòu)復(fù)雜的方程其求解結(jié)果的誤差要遠(yuǎn)比此大得多, 由于還存在著局部截?cái)嗾`差和整體截?cái)嗾`差, 因此可采取一定措施來抑制減少誤差, 盡量使結(jié)果精確.（2）使用改進(jìn)歐拉算法進(jìn)行一般求解經(jīng)過改進(jìn)歐拉算法程序計(jì)算得出在各點(diǎn)處的近似數(shù)值解及各自同精確解的誤差, 如下表3.2表3.2各分點(diǎn)(數(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)過改進(jìn)的歐拉算法所存在的誤差已大為減少, 誤差的減少主要是由于先利用了歐拉公式對(duì)的值進(jìn)行了預(yù)估, 然后又利用梯形公式對(duì)預(yù)估值作了校正, 從而在預(yù)估校正的過程中減少了誤差. 此算法有一定的優(yōu)點(diǎn), 在一些實(shí)際而且簡(jiǎn)單的工程計(jì)算中可以直接單獨(dú)應(yīng)用.（3）使用經(jīng)典RungeKutta法進(jìn)行一般求解經(jīng)過經(jīng)典RungeKutta法程序計(jì)算得出在各點(diǎn)處的近似數(shù)值解及各自同精確解的誤差, 如下表3.3表3.3各分點(diǎn)(數(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記錄的程序運(yùn)行結(jié)果來看, 所計(jì)算出來的常微分方程的數(shù)值解是非常精確的, 用浮點(diǎn)數(shù)進(jìn)行

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

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

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