常微分方程初值問題的runge-kutta解法【信息科學(xué)與技術(shù)專業(yè)】【畢業(yè)設(shè)計+文獻(xiàn)綜述+開題報告】

1、<p>　　本科畢業(yè)論文（設(shè)計）</p><p><b>　?。ā?0 　屆）</b></p><p>　　常微分方程初值問題的Runge-Kutta解法</p><p>　　所在學(xué)院 </p><p>　　專業(yè)班級 信息與計算科學(xué)

2、 </p><p>　　學(xué)生姓名 學(xué)號 </p><p>　　指導(dǎo)教師 職稱 </p><p>　　完成日期 年 月 </p><p>　　摘要：常微分方程是指只有一個自變量的微分方程。常微分方程在很多學(xué)

3、科領(lǐng)域內(nèi)有著重要的作用，如自動控制、各種電子學(xué)裝置的設(shè)計等。本文首先介紹了常微分方程以及初值問題的除Runge-Kutta解法外的數(shù)值解一些理論；接著，我們對常微分方程初值問題提出Runge-Kutta解法的構(gòu)造公式；然后，我們再對這個公式進(jìn)行誤差分析，包括相容性、穩(wěn)定性和收斂性分析；最后，通過一些數(shù)值實例，我們驗證了Runge-Kutta法的優(yōu)缺點。</p><p>　　關(guān)鍵詞：常微分方程；Runge-Kutt

4、a解法；誤差分析 </p><p>　　The Runge-Kutta Method for Initial Value Problems </p><p>　　of Ordinary Differential Equations </p><p>　　Abstract: The ordinary differential equations have only

5、one independent variable in the equation. And the ordinary differential equations play an important role in many areas, such as automatic control, electronics equipment design and so on. In this paper, we firstly introdu

6、ce some theories of the initial value problems for ordinary differential equations, except the Runge-Kutta method; Next, we construct some numerical formulations of the Runge-Kutta Method to solve the ??initial value pro

7、ble</p><p>　　Keywords: ordinary differential equations; Runge-Kutta method; error estimates</p><p><b>　　目 錄</b></p><p>　　1 緒 論1</p><p>　　1.1 問題的背景1&l

8、t;/p><p>　　1.2 常微分方程初值問題簡介1</p><p>　　2 常微分方程初值問題的一些數(shù)值解法3</p><p>　　2.1 單步法中的歐拉解法3</p><p>　　2.2 線性多步法3</p><p>　　3 Runge-Kutta解法4</p><p&

9、gt;　　3.1 Runge-Kutta解法的構(gòu)造4</p><p>　　3.2 誤差分析8</p><p>　　3.2.1 相容性8</p><p>　　3.2.2 收斂性8</p><p>　　3.2.3 穩(wěn)定性10</p><p>　　4 數(shù)值算例14</p><p

10、>　　5 結(jié) 論16</p><p><b>　　致 謝17</b></p><p><b>　　參考文獻(xiàn)18</b></p><p><b>　　附錄19</b></p><p><b>　　1 緒 論</b></p>

11、<p>　　1.1 問題的背景</p><p>　　常微分方程在微積分概念出現(xiàn)后即已出現(xiàn)。從萊布尼茲專門研究用變量變換解決一階微分方程的求解問題的“求通解”時代，到1841年劉維爾的里卡蒂方程不存在一般初等解和柯西的初值問題的提出，常微分方程轉(zhuǎn)向“求定解”時代。再到19世紀(jì)末天體力學(xué)的太陽系穩(wěn)定性問題眼界需要常微分方程解的大范圍性態(tài)，從而發(fā)展到了“求所有解問題”。最后到了20世紀(jì)六七十年代，常微分方

12、程因為計算機(jī)的發(fā)展進(jìn)入“求特殊解”階段，發(fā)現(xiàn)了具有新性質(zhì)的特殊的解和方程，如混沌（解）、奇異吸引子及孤立子等等[1]。</p><p>　　常微分方程是指只有一個自變量的微分方程。而且常微分方程在很多學(xué)科領(lǐng)域內(nèi)有著重要的作用，如自動控制、各種電子學(xué)裝置的設(shè)計、彈道的計算、飛機(jī)和導(dǎo)彈飛行的穩(wěn)定性的研究、化學(xué)反應(yīng)過程穩(wěn)定性的研究等等。常微分方程初值問題的數(shù)值求解是求近似解的一種經(jīng)典方法。</p><

13、;p>　　1.2 常微分方程初值問題簡介[2]-[9]</p><p>　　大多數(shù)學(xué)科領(lǐng)域內(nèi)的問題需要求解一個初值問題，即求解滿足給定初值條件的微分方程。我們更多的是使用逼近原問題的解的方法來逼近原方程的解。因逼近方法給出更精確的結(jié)果和實際的誤差信息，所以更經(jīng)常被使用。現(xiàn)代科學(xué)、技術(shù)、工程中的大量數(shù)學(xué)模型都可以用常微分方程的初值問題來描述，所以研究常微分方程初值問題的一些性質(zhì)、解法很有必要。一些典型的常微

14、分初值問題的數(shù)值求解問題的方法有：單步法和線性多步法。</p><p>　　在求解區(qū)間上取定節(jié)點</p><p>　　令，稱為積分網(wǎng)格的步長。用表示初值問題精確解在節(jié)點上函數(shù)值的近似值。對給定的數(shù)值積分方法，中的各個市按某一個遞推算法確定的。一個遞推算法，如果在用它計算時只用到已經(jīng)求出的諸值中的，而無須使用其余值中的任何一個，則稱此算法為單步方法。相反的則是多步法。</p>

15、<p>　　對于一階的常微分方程初值問題</p><p><b>　　， </b></p><p><b>　?。?）</b></p><p>　　的數(shù)值方法，其中是和的已知函數(shù)，為給定的初值。為使上面的解存在、唯一且連續(xù)依賴初值，則必須關(guān)于滿足Lipschitz條件：即存在常數(shù)L，使得對于任意，均成立不

16、等式</p><p>　　單步顯式公式的一般形式為</p><p><b>　　，</b></p><p><b>　　（2）</b></p><p>　　稱為增量函數(shù)。一般來說微分方程問題（1）的精確解不滿足式（2），即一般地有</p><p><b>　　定義1

17、 稱</b></p><p>　　為單步顯式公式（2）在處的截斷誤差。</p><p>　　一個求解公式的局部截斷誤差刻畫了其逼近微分方程的準(zhǔn)確程度。根據(jù)上述定義可以直接求得各式的局部截斷誤差。</p><p>　　定義2 如果一個求解公式的局部截斷誤差為，則稱該求解公式是階的，或稱具有階精度。</p><p>　　在多步

18、法中也有幾個基本的理論問題的定義，定理。</p><p>　　定義3 對于初值問題（1）的多步法</p><p>　　， （3）</p><p>　　說是相容的，如果它至少是一階的。</p><p>　　定理1 為使步法（3）相容，必需且只需</p><p>　

19、　2 常微分方程初值問題的一些數(shù)值解法</p><p>　　下面將先簡單的介紹一下除Runge-Kutta法外的一些數(shù)值方法。</p><p>　　2.1 單步法中的歐拉解法[4]-[9]</p><p>　　單步法中最簡單的數(shù)值解法是歐拉法。歐拉法的基本思想是由推導(dǎo)出下一個節(jié)點，其推導(dǎo)公式為</p><p><b>　　，.&

20、lt;/b></p><p>　　其局部截斷誤差為。歐拉法的幾何意義是用一條過的折線來近似替代過的積分曲線。</p><p>　　雖然這個方法的計算很簡單，但精確度較低，所以有了精確度較高的梯形方法。推導(dǎo)公式為，，它的局部截斷誤差為。此為隱型格式，迭代求解所需計算量大。</p><p>　　結(jié)合上述兩種方法的優(yōu)點，從而得到改進(jìn)歐拉方法。其計算公式為， 。

21、</p><p>　　2.2 線性多步法[8][10]</p><p>　　常用的多步法主要有Adams法和Milne法，其中Adams法又包括顯式Adams法和隱式Adams法。顯式Adams法的線性步公式為</p><p><b>　　其中</b></p><p><b>　　， 。</b>

22、</p><p>　　且隱式Adams公式為</p><p><b>　　其中</b></p><p><b>　　， </b></p><p>　　3 Runge-Kutta解法</p><p>　　3.1 Runge-Kutta解法的構(gòu)造[5]-[15]</

23、p><p>　　用Taylor級數(shù)法的公式，可以構(gòu)造高階單步法，但增量函數(shù)是用的各階導(dǎo)數(shù)表示的，通常不易計算。由于函數(shù)在一點的導(dǎo)數(shù)值可以用該點附近若干點的函數(shù)值近似表示，因此可以把Taylor級數(shù)法中的增量函數(shù)改為在一些點上函數(shù)值的組合，然后用Taylor展開確定待定的系數(shù)，使方法達(dá)到一定的階。這就是Runge-Kutta法的基本思想?；蛘呤窃诿總€區(qū)間多取幾點，將它們的斜率加權(quán)平均作為平均斜率，就有可能構(gòu)造出精度更高

24、的計算公式。通過構(gòu)造高階單步法來提高精度，而較高的精度意味著計算結(jié)果更加精確，誤差會隨著的減小而迅速減小。</p><p><b>　　現(xiàn)在假設(shè)</b></p><p><b>　　（4）</b></p><p><b>　　其中</b></p><p><b>　　

25、，</b></p><p>　　， （5）</p><p>　　此類遞推算法為顯示N級Runge-Kutta方法。上述算法公式中的系數(shù)，和希望如此確定，使得（4）式的Taylor展開式中所有冪次不超過的那些項的系數(shù)與Taylor展開式 （6） （其中，）中的相應(yīng)項的系數(shù)相等。</p&g

26、t;<p>　　當(dāng)N=1時，則可知=1，此時式（4）為歐拉公式。</p><p>　　當(dāng)N=2時，確定系數(shù)，和使得與兩者的Taylor展開中含項的系數(shù)，，對盡可能大的相等。</p><p>　　利用二元的Taylor展開公式</p><p><b>　　可將（4）式表示為</b></p><p><b

27、>　?。?）</b></p><p>　　再將Taylor展式（4）改成</p><p><b>　　（8）</b></p><p>　　比較（7）和（8）式中項的系數(shù)，得到</p><p>　　， ， （9）</p><

28、p>　　因為有4個待定系數(shù)，可是就只有3個方程數(shù)目。因此具有無窮多個數(shù)組滿足條件（9），即存在無窮多個二級Runge-Kutta方法，它們的截斷誤差。</p><p><b>　　若取，便得到</b></p><p><b>　　或</b></p><p><b>　　這是改進(jìn)歐拉公式。</b>

29、</p><p>　　因為二級的RK方法最多是二階的，若想得到更高階的方法，必須增加級數(shù)，即取較大的N。當(dāng)N較大時，計算也就會更加的復(fù)雜。</p><p>　　當(dāng)N=3時，和N=2的推導(dǎo)過程相相似，當(dāng)時，得到</p><p><b>　　或當(dāng)時，有</b></p><p>　　同理，N=4時，構(gòu)造經(jīng)典的四階RK方法，其計

30、算公式為</p><p><b>　?。?0）</b></p><p><b>　　其局部截斷誤差為。</b></p><p>　　但是，以上的方法還是不夠的。在涉及初值問題數(shù)值解的實際問題中，總要先設(shè)定一個精度允許值，數(shù)值解與真解的偏離不能超過這個允許值。因此，步長還受到這個允許值的控制，常常會發(fā)生這樣的情形：僅在求解的

31、某一部分需要小的步長，而在其余部分則用較大的步長就可以了。</p><p>　　所以，各種用于算法中自動地調(diào)整步長的方法（稱為自適應(yīng)方法）就得到了發(fā)展。就如在四階的Runge-Kutta方法中，為了從到，我們可以用步長為做一步，也可以用步長為做兩步，如果沒有截斷誤差，它們求出的的值應(yīng)是相同的。但事實上它們不會相同，所以其差異可看成是對局部截斷誤差的一個估計：如果這個差異不超過預(yù)先給定的允許值，那么步長就是滿意的；

32、如果它超過允許值，就將減半；如果它遠(yuǎn)小于允許值，就將步長加倍。</p><p>　　雖然這樣的過程很容易計算，但要浪費(fèi)大量的時間，因而不宜提倡。</p><p>　　所以給出一個Runge-Kutta方法的漂亮的自適應(yīng)格式:</p><p><b>　?。?1）</b></p><p>　　它要求計算五個函數(shù)值，但只是一

33、個四階方法，因此單獨(dú)使用并無好處，但它可以通過增加一個函數(shù)值（而不需要改變已有的五個函數(shù)值）來得到一個更高一階的方法：</p><p>　　同（11）式 （12）</p><p>　　從四階與五階過程中得到的的值之間的差，是對四階方法中局部截斷誤差的一個估計。換句話說，通過求六個函數(shù)值和簡單的四則運(yùn)算（系數(shù)可以預(yù)先計算

34、好存放起來），我們得到了一個四階的近似值連同一個誤差估計！計算機(jī)只要比較這個近似值和允許誤差界，就可以自動地決定在下一步計算中，到底是保持還是縮放當(dāng)前正在使用的步長。</p><p>　　由上可知，由于RK方法是通過對函數(shù)的Taylor展開來構(gòu)造的，因此高階RK方法要求微分方程的解具有高階導(dǎo)數(shù)。如果解的光滑性較差，即使用高階方法也不可能得到高精度的數(shù)值解，此時 應(yīng)當(dāng)用低階方法。或者說從這些上述公式的討論中可知，級

35、數(shù)N表示計算函數(shù)值的次數(shù)，當(dāng)時，N級的RK方法的階數(shù)，即計算函數(shù)值的次數(shù)越多，RK方法的精度即階數(shù)越高。但是當(dāng)時，情況就大不一樣了，J.C.Butcher指出：當(dāng)時，；當(dāng)時，；而且，。由此可見四級以上的RK方法，計算函數(shù)值的工作量增加較快，而階數(shù)即精度提高較慢，故而，在實際問題中比較常用的是四級RK方法。</p><p>　　3.2 誤差分析[9]-[15]</p><p>　　3.2.

36、1 相容性</p><p>　　定義4 當(dāng)，時，稱數(shù)值方法是相容的。</p><p><b>　　由于</b></p><p><b>　　于是</b></p><p>　　故有結(jié)論：R-K法相容。即p階的Runge-Kutta方法是相容的。</p><p>　　3.2.

37、2 收斂性</p><p>　　由于截斷誤差的存在，計數(shù)值只是精確解的一個近似。故我們可以相當(dāng)于有這樣的問題：差是不是隨著步長的減小而減小，這就是收斂性問題。</p><p>　　定義5 當(dāng)時，，，稱數(shù)值方法是收斂的。</p><p><b>　　令 ，則</b></p><p><b>　　-） <

38、;/b></p><p><b>　　于是得到結(jié)論：</b></p><p>　　當(dāng)時，，且有誤差估計</p><p><b>　　。</b></p><p>　　即也有由于，當(dāng)時，，當(dāng)步長減小時，從計算到點的步數(shù)就增大，從而推出下面的定理。</p><p>　　定理2

39、 若公式（2）是階方法，且增量函數(shù)關(guān)于滿足Lipschitz條件，則該方法收斂。</p><p>　　對于Runge-Kutta方法，當(dāng)函數(shù)關(guān)于滿足Lipschitz條件，可用數(shù)學(xué)歸納法證明所有的關(guān)于滿足Lipschitz條件，從而增量函數(shù)關(guān)于滿足Lipschitz條件，所以Runge-Kutta方法收斂。而且對于這類滿足Lipschitz條件的單步法而言，相容條件是收斂條件的充分必要條件。</p>

40、<p>　　3.2.3 穩(wěn)定性</p><p>　　下面討論穩(wěn)定性問題，更準(zhǔn)確地說就是數(shù)值解關(guān)于初值的連續(xù)依賴性問題。</p><p>　　定義6 如果一個單步法滿足：存在正常數(shù)和，使得用任意步長以及任意兩個初值和計算得到的兩組數(shù)值和成立</p><p><b>　　，</b></p><p>　　則稱

41、該方法是穩(wěn)定的，其中。</p><p>　　定理3 若增量函數(shù)關(guān)于滿足Lipschitz條件，則公式（2）是穩(wěn)定的。</p><p><b>　　令 </b></p><p><b>　　-） </b></p><p><b>　　于是得到</b></p>

42、<p>　　其中是的次多項式。一般地，</p><p><b>　　， </b></p><p><b>　　于是</b></p><p><b>　　記，則</b></p><p>　　即Runge-Kutta方法穩(wěn)定。</p><p>　　

43、以上是在的條件下才滿足的。故可知，歐拉法，改進(jìn)歐拉法和Runge-Kutta方法都是穩(wěn)定的。然而，第一，我們實際上是取有限的固定步長進(jìn)行計算的，它并不能隨意縮小至與零充分接近。第二，即使可以任意縮小，我們也不會那么做，這將導(dǎo)致計算工作亮的急劇增加。第三，更重要的是，由于舍入誤差在計算過程中的產(chǎn)生、傳播和放大，小的步長可能導(dǎo)致計算結(jié)果嚴(yán)重失真甚至面目全非。因此在現(xiàn)實中我們關(guān)心的恰恰是對一個方法來講，它允許的最大步長是多少。顯然，這個最大步

44、長與所給方程的性質(zhì)有關(guān)。所以就又有了使任何一步產(chǎn)生的誤差在以后過程中足夠減少的絕對穩(wěn)定性。</p><p>　　定義7 對于Runge-Kutta方法的</p><p><b>　　當(dāng)即充分小時，</b></p><p>　　Runge-Kutta方法穩(wěn)定滿足根條件</p><p><b>　　當(dāng)不趨近0時，

45、</b></p><p><b>　　-） </b></p><p><b>　　一般地，記 ，，則</b></p><p>　　其中是的次齊次多項式，且其系數(shù)之和為。于是</p><p><b>　　即</b></p><p>　　Run

46、ge-Kutta方法絕對穩(wěn)定滿足嚴(yán)格根條件。即求，使</p><p><b>　　滿足：</b></p><p>　　這時，所求的的變動范圍就是絕對穩(wěn)定區(qū)域。</p><p>　　下表3-1給出了高階單步方法在實數(shù)域中的絕對穩(wěn)定區(qū)間：</p><p>　　表3-1絕對穩(wěn)定區(qū)間</p><p>　　

47、由上表可以發(fā)現(xiàn)，當(dāng)增加時它的絕對穩(wěn)定區(qū)域還有所增大[15]。</p><p>　　故由此可知，當(dāng)滿足時，N級階的Runge-Kutta方法(=1,2,3,4)是絕對穩(wěn)定的。</p><p>　　如果一個方法的絕對穩(wěn)定區(qū)域比較大，則可說它的絕對穩(wěn)定性比較好，且隱式方法的絕對穩(wěn)定性比顯式方法好；對于Runge-Kutta方法，階愈高絕對穩(wěn)定性愈好。</p><p>&l

48、t;b>　　4 數(shù)值算例</b></p><p>　　因為現(xiàn)在無法比較出4階的Runge-Kutta方法與其他階數(shù)求解方法的優(yōu)劣，所以這里引進(jìn)比較低階Runge-Kutta方法的一個度量描述如下：</p><p>　　4階Runge-Kutta方法每步需要4次求值，所以如果它有更優(yōu)的話，它應(yīng)能比1/4步長的歐拉法給出更精確地答案。類似地，如果4階Runge-Kutta方

49、法比二階Runge-Kutta方法優(yōu)越的話，則步長為的這種方法應(yīng)該比步長為的二階方法給出更精確地結(jié)果，因為4階方法每步需要二倍次數(shù)的求值。</p><p>　　算例4-1 對于問題</p><p><b>　　，，，</b></p><p>　　用四階的Runge-Kutta方法解上式問題。</p><p>　　解：對

50、于上式我們采用不同階數(shù)的Runge-Kutta方法，即，，，時，令步長，從而得出各種不同階數(shù)的Runge-Kutta方法的相對誤差，如下表4-1所示。</p><p><b>　　表4-1 </b></p><p>　　由以上的表可以觀察到4階Runge-Kutta方法的誤差是最小的。即比其它的階數(shù)Runge-Kutta方法有更高的精確度。</p>&l

51、t;p>　　算例4-2 對于問題</p><p><b>　　，，</b></p><p>　　分別取步長為（i=1,2），求解。</p><p>　　解：當(dāng)時，不屬于絕對穩(wěn)定區(qū)域，而步長時，相應(yīng)的屬于絕對穩(wěn)定區(qū)域</p><p>　　計算結(jié)果如下表4-2所示。</p><p><b

52、>　　表4-2</b></p><p>　　從上表中，我們可以觀察到當(dāng)時，其相對誤差是呈現(xiàn)出收斂與常數(shù)的趨勢。但當(dāng)時，其相對誤差就呈現(xiàn)發(fā)散的趨勢，顯然這是由于此時的不屬于絕對穩(wěn)定區(qū)域所引起的，故我們可以知道絕對穩(wěn)定性的重要。</p><p><b>　　5 結(jié) 論</b></p><p>　　本文先介紹了常微分方程的初值問題

53、的一些數(shù)值解法，在詳細(xì)介紹了單步法中的Runge-Kutta解法。常微分方程是指只有一個自變量的微分方程，而且常微分方程在很多學(xué)科領(lǐng)域內(nèi)有著重要的作用。常微分方程初值問題的數(shù)值求解是求近似解的一種經(jīng)典方法，Runge-Kutta方法是目前工程上求常微分方程數(shù)值解的基本方法。它的基本思想是在每個區(qū)間多取幾點，將它們的斜率加權(quán)平均作為平均斜率，就有可能構(gòu)造出精度更高的計算公式。</p><p>　　通過對Runge-

54、Kutta方法的構(gòu)造，通過對Runge-Kutta方法的研究，得出了Runge-Kutta法的相容性、收斂性和穩(wěn)定性分析，通過程序（MATLAB語言）編程，給出數(shù)值例子，尤其是4階的Runge-Kutta公式，驗證了理論分析。常微分方程初值問題的數(shù)值解在數(shù)學(xué)意義上的存在性可以在一般的條件下得到證明，而從實際應(yīng)用上來講，人們需要的往往并不是解在數(shù)學(xué)上的存在性，而是在我們所關(guān)心的某個定義范圍內(nèi)，對應(yīng)于某些特定的自變量的取值，尋求數(shù)值解，而這

55、也是為什么要對其數(shù)值計算法中Runge-Kutta法的研究。</p><p><b>　　參考文獻(xiàn)</b></p><p>　　[1] 王高雄，朱思銘，王壽松．常微分方程[M]．北京：高等教育出版社，2006．</p><p>　　[2] [美]Richard L.Burden，J.Douglas Faires． 數(shù)值分析[[M]．北京：高等教

56、育出版社，2005．</p><p>　　[3] 中國大百科全書·數(shù)學(xué)．北京：中國大百科全書出版社，1988．</p><p>　　[4] 余徳浩，湯華中．微分方程數(shù)值解法[M]．北京：科學(xué)出版社，2003．</p><p>　　[5] 胡健偉，湯懷民．微分方程數(shù)值解法[M]．北京：科學(xué)出版社，1999．</p><p>　　[6]

57、 李榮華，馮果忱．微分方程數(shù)值解法[M]．北京：高等教育出版社，1997．</p><p>　　[7] Arieh Iserles．微分方程數(shù)值分析基礎(chǔ)教程[M]．北京：清華大學(xué)出版社，2005．</p><p>　　[8] 李忠杰．常微分方程初值問題RK法和多步法[J]．黑龍江科技信息，科教文化，（1）：199．</p><p>　　[9] 李瑞遐，何志慶．微分方

58、程數(shù)值方法[M]．上海：華東理工大學(xué)出版社，2005．</p><p>　　[10] 孫維夫，姜云義．常微分方程初值問題RK法和多步法[J]．煙臺職業(yè)學(xué)院學(xué)報，2010，（2）：84-88．</p><p>　　[11] 劉利斌，王偉．四階Runge-Kutta算法的優(yōu)化分析[J]．成都大學(xué)學(xué)報（自然科學(xué)版），2007，（1）：19-21．</p><p>　　[1

59、2] 林群．微分方程數(shù)值解法基礎(chǔ)教程[M]．北京：科學(xué)出版社，2003．</p><p>　　[13] 馬明書．龍格-庫塔法的級數(shù)和階數(shù)[J]．河南師范大學(xué)河南機(jī)專學(xué)報，1994，（2）：5-7．</p><p>　　[14] 吳怡，張向輝，熊超西．三階Runge-Kutta最優(yōu)算法[J]．河北秦皇島職業(yè)技術(shù)學(xué)院基礎(chǔ)部，2005，（18）：57-61．</p><p&g

60、t;　　[15] 李立康，於崇華．微分方程數(shù)值方法[M]．上海：復(fù)旦大學(xué)出版社，2003．</p><p><b>　　附錄</b></p><p>　　附錄1：例題4.1的編程</p><p><b>　　函數(shù)的表達(dá)式：</b></p><p>　　function f=fun01(x,y)<

61、;/p><p>　　f=y-2*x/y;</p><p>　　function ex01()</p><p>　　x0=0;xt=1;</p><p><b>　　y0=1;</b></p><p>　　PointNum=16;</p><p>　　[x1,y1] = MyEu

62、ler1D('fun01',x0,xt,y0,PointNum+1);</p><p>　　[x2,y2] = MyRunge_Kutta2('fun01',x0,xt,y0,PointNum);</p><p>　　[x3,y3] = MyRunge_Kutta3('fun01',x0,xt,y0,PointNum);</p>

63、<p>　　[x4,y4] = MyRunge_Kutta4('fun01',x0,xt,y0,PointNum);</p><p>　　y_exact=(2*x3+1).^(1/2);</p><p>　　Err1=abs(y1-y_exact)./abs(y_exact)</p><p>　　Err2=abs(y2-y_exact).

64、/abs(y_exact)</p><p>　　Err3=abs(y3-y_exact)./abs(y_exact)</p><p>　　Err4=abs(y4-y_exact)./abs(y_exact)</p><p>　　function [outx,outy]=MyEuler1D(fun,x0,xt,y0,PointNum)</p><p

65、>　　%MyEuler1D 用歐拉方法解微分方程</p><p>　　%fun 表示f（x，y）</p><p>　　%x0,xt表示自變量的初值和終值</p><p>　　%y0表示函數(shù)在x0處的值,可以是向量</p><p>　　%PointNum表示自變量在[x0,xt]上取的點數(shù),包括端點.</p><p&g

66、t;　　if nargin<5 | PointNum<=0 %如果函數(shù)僅輸入4個參數(shù)值，則PointNum默認(rèn)值為101</p><p>　　PointNum=101;</p><p><b>　　end</b></p><p>　　if nargin<4 %y0默認(rèn)值為0</p><p><b&

67、gt;　　y0=0;</b></p><p><b>　　end</b></p><p>　　h=(xt-x0)/(PointNum-1);%計算步長h</p><p>　　%x=x0:h:xt;</p><p>　　x=x0+[0:PointNum-1]'*h;%自變量數(shù)組</p>&

68、lt;p>　　y(1,:)= y0;%將輸入存為行向量，輸入為列向量形式</p><p><b>　　%pause</b></p><p>　　for k = 1:PointNum-1</p><p>　　f=feval(fun,x(k),y(k,:));%計算f(x,y)在每個迭代點的值</p><p>　　%

69、 f=f(:)',</p><p>　　y(k + 1,:) =y(k,:) +h*f; %對于所取的點x迭代計算y值</p><p><b>　　end</b></p><p><b>　　outy=y;</b></p><p><b>　　outx=x;</b>

70、</p><p>　　%plot(x,y)%畫出方程解的函數(shù)圖</p><p>　　function [x,y] = MyRunge_Kutta2(fun,x0,xt,y0,PointNum,varargin)</p><p>　　%二階Runge-Kutta 方法解微分方程形為 y’(t) = f(x,y(x))</p><p>　　%此程

71、序可解高階的微分方程。只要將其形式寫為上述微分方程的向量形式</p><p>　　% x范圍為[x0,xt],初值為 y0, PointNum為離散點數(shù),varargin為可選輸入項可傳適當(dāng)參數(shù)給函數(shù)f(x,y)</p><p>　　if nargin < 4 | PointNum <= 0</p><p>　　PointNum= 100; </p

72、><p><b>　　end</b></p><p>　　if nargin < 3</p><p><b>　　y0 = 0; </b></p><p><b>　　end</b></p><p>　　y(1,:) = y0(:)'; %初值

73、存為行向量形式</p><p>　　h = (xt-x0)/(PointNum-1); %計算步長</p><p>　　x = x0+[0:PointNum]'*h; %得x向量值</p><p>　　for k = 1:PointNum %迭代計算</p><p>　　f1 = h*feval(fun,x(k),y(k,:),va

74、rargin{:}); </p><p>　　f1 = f1(:)'; %得公式中k1</p><p>　　f2 = h*feval(fun,x(k) + h/2,y(k,:) + f1/2,varargin{:}); </p><p>　　f2 = f2(:)'; %得公式中k2</p><p>　　y(k + 1,:) =

75、 y(k,:) + f2; %得y(n+1)</p><p><b>　　end</b></p><p>　　function [x,y] = MyRunge_Kutta3(fun,x0,xt,y0,PointNum,varargin)</p><p>　　%三階Runge-Kutta 方法解微分方程形為 y’(t) = f(x,y(x))&l

76、t;/p><p>　　%此程序可解高階的微分方程。只要將其形式寫為上述微分方程的向量形式</p><p>　　% x范圍為[x0,xt],初值為 y0, PointNum為離散點數(shù),varargin為可選輸入項可傳適當(dāng)參數(shù)給函數(shù)f(x,y)</p><p>　　if nargin < 4 | PointNum <= 0</p><p>

77、;　　PointNum= 100; </p><p><b>　　end</b></p><p>　　if nargin < 3</p><p><b>　　y0 = 0; </b></p><p><b>　　end</b></p><p>　　

78、y(1,:) = y0(:)'; %初值存為行向量形式</p><p>　　h = (xt-x0)/(PointNum-1); %計算步長</p><p>　　x = x0+[0:PointNum]'*h; %得x向量值</p><p>　　for k = 1:PointNum %迭代計算</p><p>　　f1 = h*

79、feval(fun,x(k),y(k,:),varargin{:}); </p><p>　　f1 = f1(:)'; %得公式中k1</p><p>　　f2 = h*feval(fun,x(k) + h/3,y(k,:) + f1/3,varargin{:}); </p><p>　　f2 = f2(:)'; %得公式中k2</p>

80、<p>　　f3 = h*feval(fun,x(k) + 2*h/3,y(k,:) + 2*f2/3,varargin{:}); </p><p>　　f3 = f3(:)'; %得公式中k3</p><p>　　y(k + 1,:) = y(k,:) + (f1 + 3*f3)/4; %得y(n+1)</p><p><b>　　e

81、nd</b></p><p>　　function [x,y] = MyRunge_Kutta4(fun,x0,xt,y0,PointNum,varargin)</p><p>　　%四階Runge-Kutta 方法解微分方程形為 y’(t) = f(x,y(x))</p><p>　　%此程序可解高階的微分方程。只要將其形式寫為上述微分方程的向量形式&

82、lt;/p><p>　　% x范圍為[x0,xt],初值為 y0, PointNum為離散點數(shù),varargin為可選輸入項可傳適當(dāng)參數(shù)給函數(shù)f(x,y)</p><p>　　if nargin < 4 | PointNum <= 0</p><p>　　PointNum= 100; </p><p><b>　　end&l

83、t;/b></p><p>　　if nargin < 3</p><p><b>　　y0 = 0; </b></p><p><b>　　end</b></p><p>　　y(1,:) = y0(:)'; %初值存為行向量形式</p><p>　　h

84、 = (xt-x0)/(PointNum-1); %計算步長</p><p>　　x = x0+[0:PointNum]'*h; %得x向量值</p><p>　　for k = 1:PointNum %迭代計算</p><p>　　f1 = h*feval(fun,x(k),y(k,:),varargin{:}); </p><p&g

85、t;　　f1 = f1(:)'; %得公式中k1</p><p>　　f2 = h*feval(fun,x(k) + h/2,y(k,:) + f1/2,varargin{:}); </p><p>　　f2 = f2(:)'; %得公式中k2</p><p>　　f3 = h*feval(fun,x(k) + h/2,y(k,:) + f2/2,v

86、arargin{:}); </p><p>　　f3 = f3(:)'; %得公式中k3</p><p>　　f4 = h*feval(fun,x(k) + h,y(k,:) + f3,varargin{:}); </p><p>　　f4 = f4(:)'; %得公式中k4</p><p>　　y(k + 1,:) = y(

87、k,:) + (f1 + 2*(f2 + f3) + f4)/6; %得y(n+1)</p><p><b>　　end</b></p><p>　　附錄2：例題4.2的編程</p><p>　　function f=fun02(x,y)</p><p><b>　　f=-24*y;</b><

88、/p><p>　　function ex02()</p><p>　　x0=0;xt=1;</p><p><b>　　y0=1;</b></p><p>　　PointNum=8;</p><p>　　[x,y] = MyRunge_Kutta4('fun02',x0,xt,y0,P

89、ointNum);</p><p>　　x=x0+[0:PointNum]'*1/8;</p><p>　　y_exact=exp(-24*x);</p><p>　　plot(x,y,'o-',x,y_exact,'+:')</p><p>　　Err=abs(y-y_exact)./abs(y_ex

90、act);</p><p><b>　　Err</b></p><p>　　function ex02()</p><p>　　x0=0;xt=1;</p><p><b>　　y0=1;</b></p><p>　　PointNum=32;</p><p

91、>　　[x,y] = MyRunge_Kutta4('fun02',x0,xt,y0,PointNum);</p><p>　　x=x0+[0:PointNum]'*1/32;</p><p>　　y_exact=exp(-24*x);</p><p>　　plot(x,y,'o-',x,y_exact,'+:&

92、#39;)</p><p><b>　　hold on;</b></p><p>　　Err=abs(y-y_exact)./abs(y_exact);</p><p><b>　　Err</b></p><p>　　err=Err(4:4:32)</p><p><b

93、>　　文獻(xiàn)綜述</b></p><p>　　常微分方程初值問題的Runge-Kutta解法 　　　　　　　　　　　　　　　　　</p><p><b>　　一、前言部分</b></p><p>　　常微分方程在很多學(xué)科領(lǐng)域內(nèi)有著重要的作用，如自動控制、各種電子學(xué)裝置的設(shè)計、彈道的計算、飛機(jī)和導(dǎo)彈飛行的穩(wěn)定性的研究、化學(xué)反應(yīng)

94、過程穩(wěn)定性的研究等等，這些問題都可以化為求微分方程的解，或者化為研究解的性質(zhì)的問題。這些問題都包含某個變量關(guān)于另一個變量的變化。大多數(shù)這樣的問題需要求解一個初值問題，即求解滿足給定初值條件的微分方程。我們更多的是使用逼近原問題的解的方法來逼近原方程的解。因逼近方法給出更精確的結(jié)果和實際的誤差信息，所以更經(jīng)常被使用。一些典型的常微分初值問題的數(shù)值求解問題的方法有：單步法和線性多步法。</p><p>　　在求解區(qū)間

95、上取定節(jié)點</p><p>　　令，稱為積分網(wǎng)格的步長。用表示初值問題精確解在節(jié)點上函數(shù)值的近似值。對給定的數(shù)值積分方法，中的各個市按某一個遞推算法確定的。一個遞推算法，如果在用它計算時只用到已經(jīng)求出的諸值中的，而無須使用其余值中的任何一個，則稱此算法為單步方法。相反的則是多步法。單步法主要有歐拉法和Runge-Kutta法，多步法主要有Adams法和Milne法等。</p><p>　　

96、本文綜述常微分初值問題初值問題的數(shù)值解法及其誤差估計（相容性、穩(wěn)定性和收斂性分析），重點介紹了Runge-Kutta法。</p><p><b>　　二、主題部分</b></p><p>　　2.1 常微分方程的初值問題概述</p><p>　　常微分方程在微積分概念出現(xiàn)后即已出現(xiàn)。從萊布尼茲專門研究用變量變換解決一階微分方程的求解問題的“求通

97、解”時代，到1841年劉維爾的里卡蒂方程不存在一般初等解和柯西的初值問題的提出，常微分方程轉(zhuǎn)向“求定解”時代。再到19世紀(jì)末天體力學(xué)的太陽系穩(wěn)定性問題眼界需要常微分方程解的大范圍性態(tài)，從而發(fā)展到了“求所有解問題”。最后到了20世紀(jì)六七十年代，常微分方程因為計算機(jī)的發(fā)展進(jìn)入“求特殊解”階段，發(fā)現(xiàn)了具有新性質(zhì)的特殊的解和方程，如混沌（解）、奇異吸引子及孤立子等等。</p><p>　　在數(shù)學(xué)分析（微積分）中研究了變量

98、的各種函數(shù)及函數(shù)的微分與積分。如函數(shù)未知，但知道變量與函數(shù)的代數(shù)關(guān)系式，便組成代數(shù)方程，通過求解代數(shù)方程解出未知函數(shù)。同樣，如果知道自變量、未知函數(shù)及函數(shù)的導(dǎo)數(shù)（或微分）組成的關(guān)系式，得到的便是微分方程，通過求解微分方程求出未知函數(shù)。自變量只有一個的微分方程稱為常微分方程。常微分方程是數(shù)學(xué)分析或基礎(chǔ)數(shù)學(xué)的一個組成部分，在整個數(shù)學(xué)大廈中占據(jù)著重要位置。</p><p>　　在常微分方程中，我們知道只有少數(shù)簡單的常微

99、分方程才能用初等積分法求出它們的精確解析解，多數(shù)情形只能用近似解法求其近似解。近似解法可以分成兩類：解析方法和數(shù)值方法。其中前者尋求解的近似表達(dá)式，后者則是計算微分方程解在求解區(qū)域中一些離散點上的近似值。</p><p>　　對于一階的常微分方程初值問題</p><p><b>　　， </b></p><p><b>　　（

100、1）</b></p><p>　　的數(shù)值方法，其中是和的已知函數(shù)，為給定的初值。為使上面的解存在、唯一且連續(xù)依賴初值，則必須關(guān)于滿足Lipschitz條件：即存在常數(shù)L，使得對于任意，均成立不等式</p><p>　　單步顯式公式的一般形式為</p><p><b>　　，</b></p><p><b

101、>　?。?）</b></p><p>　　稱為增量函數(shù)。一般來說微分方程問題（1）的精確解不滿足式（2），即一般地有</p><p><b>　　定義1 稱</b></p><p>　　為單步顯式公式（2）在處的截斷誤差。</p><p>　　一個求解公式的局部截斷誤差刻畫了其逼近微分方程的準(zhǔn)確程度

102、。根據(jù)上述定義可以直接求的各式的局部截斷誤差。例如歐拉公式的局部截斷誤差為</p><p>　　定義2 如果一個求解公式的局部截斷誤差為，則稱該求解公式是階的，或稱具有階精度。</p><p>　　2.2 常微分方程的初值問題的數(shù)值方法</p><p>　　下面將先簡單的介紹一下除Runge-Kutta法外的一些數(shù)值方法。</p><p>

103、;　　2.2.1 單步法中的歐拉法</p><p>　　單步法中最簡單的數(shù)值解法是歐拉法。歐拉法的基本思想是由推導(dǎo)出下一個節(jié)點，其推導(dǎo)公式為</p><p><b>　　，.</b></p><p>　　其局部截斷誤差為。歐拉法的幾何意義是用一條過的折線來近似替代過的積分曲線。</p><p>　　雖然這個方法的計算很

104、簡單，但精確度較低，所以有了精確度較高的梯形方法。推導(dǎo)公式為，，它的局部截斷誤差為。此為隱型格式，迭代求解所需計算量大。</p><p>　　結(jié)合上述兩種方法的優(yōu)點，從而得到改進(jìn)歐拉方法。其計算公式為， 。</p><p>　　2.2.2 線性多步法</p><p>　　常用的多步法主要有Adams法和Milne法，其中Adams法又包括顯式Adams法和隱

105、式Adams法。顯式Adams法的線性步公式為</p><p><b>　　其中</b></p><p><b>　　， 。</b></p><p>　　且隱式Adams公式為</p><p><b>　　其中</b></p><p><b>

106、;　　， </b></p><p>　　2.3 Runge-Kutta法</p><p>　　2.3.1 Runge-Kutta法的構(gòu)造</p><p>　　用Taylor級數(shù)法的公式，可以構(gòu)造高階單步法，但增量函數(shù)是用的各階導(dǎo)數(shù)表示的，通常不易計算。由于函數(shù)在一點的導(dǎo)數(shù)值可以用該點附近若干點的函數(shù)值近似表示，因此可以把Taylor級數(shù)法中的增量

107、函數(shù)改為在一些點上函數(shù)值的組合，然后用Taylor展開確定待定的系數(shù)，使方法達(dá)到一定的階。這就是Runge-Kutta法的基本思想?；蛘呤窃诿總€區(qū)間多取幾點，將它們的斜率加權(quán)平均作為平均斜率，就有可能構(gòu)造出精度更高的計算公式。通過構(gòu)造高階單步法來提高精度，而較高的精度意味著計算結(jié)果更加精確，誤差會隨著的減小而迅速減小。</p><p><b>　　現(xiàn)在假設(shè)</b></p>&l

108、t;p><b>　　（3）</b></p><p><b>　　其中</b></p><p><b>　　，</b></p><p>　　， （4）</p><p>　　此類遞推算法為顯示n級Runge-Kutta方法。上述算法公式中的系數(shù)，和希望如此確定，使得

109、（3）式的Taylor展開式中所有冪次不超過的那些項的系數(shù)與Taylor展開式 （5）（其中，）中的相應(yīng)項的系數(shù)相等。</p><p>　　當(dāng)n=1時，則可知=1，此時式（3）為歐拉公式。</p><p>　　當(dāng)n=2時，確定系數(shù)，和使得與兩者的Taylor展開中含項的系數(shù)，，對盡可能大的相等。</p><p>　　利用二元的Taylor

110、展開公式</p><p><b>　　可將（3）式表示為</b></p><p><b>　?。?）</b></p><p>　　再將Taylor展式（3）改成</p><p><b>　　（7）</b></p><p>　　比較（6）和（7）式中項的系

111、數(shù)，得到</p><p>　　， ， （8）</p><p>　　因為有4個待定系數(shù)，可是就只有3個方程數(shù)目。因此具有無窮多個數(shù)組滿足條件（8），即存在無窮多個二級Runge-Kutta方法，它們的截斷誤差。</p><p><b>　　若取，便得到</b></p><p&

112、gt;<b>　　或</b></p><p><b>　　這是改進(jìn)歐拉公式。</b></p><p>　　因為二級的RK方法最多是二階的，若想得到更高階的方法，必須增加級數(shù)，即取較大的n。當(dāng)n較大時，計算也就會更加的復(fù)雜。</p><p>　　當(dāng)n=3時，和n=2的推導(dǎo)過程相相似，當(dāng)時，得到</p><

113、p><b>　　或當(dāng)時，有</b></p><p>　　同理，n=4時，構(gòu)造經(jīng)典的四階RK方法，其計算公式為</p><p><b>　?。?）</b></p><p><b>　　其局部截斷誤差為。</b></p><p>　　由于RK方法是通過對函數(shù)的Taylor展開

114、來構(gòu)造的，因此高階RK方法要求微分方程的解具有高階導(dǎo)數(shù)。如果解的光滑性較差，即使用高階方法也不可能得到高精度的數(shù)值解，此時 應(yīng)當(dāng)用低階方法?；蛘哒f從這些上述公式的討論中可知，級數(shù)n表示計算函數(shù)值的次數(shù)，當(dāng)時，n級的RK方法的階數(shù)，即計算函數(shù)值的次數(shù)越多，RK方法的精度即階數(shù)越高。但是當(dāng)時，情況就大不一樣了，J.C.Butcher指出：當(dāng)時，；當(dāng)時，；而且，。由此可見四級以上的RK方法，計算函數(shù)值的工作量增加較快，而階數(shù)即精度提高較慢，故

115、而，在實際問題中比較常用的是四級RK方法。</p><p>　　2.3.2 Runge-Kutta法的誤差分析</p><p><b>　　I、相容性</b></p><p>　　當(dāng)，時，稱數(shù)值方法是相容的。</p><p><b>　　由于</b></p><p><

116、;b>　　于是</b></p><p>　　故有結(jié)論：R-K法相容。即p階的Runge-Kutta方法是相容的。</p><p><b>　　II、收斂性</b></p><p>　　由于截斷誤差的存在，計數(shù)值只是精確解的一個近似。故我們可以相當(dāng)這樣的問題：差是不是隨著步長的減小而減小。這就是收斂性問題。</p>

117、<p>　　定義3 如果一個單步法滿足：對中任意一個固定的點成立</p><p>　　則稱該方法是收斂的。</p><p>　　由于，當(dāng)時，，也即當(dāng)步長減小時，從計算到點的步數(shù)就增大，故有下面的定理。</p><p>　　定理1 若公式（2）是階方法，且增量函數(shù)關(guān)于滿足Lipschitz條件，則該方法收斂。</p><p>　

118、　對于Runge-Kutta方法，當(dāng)函數(shù)關(guān)于滿足Lipschitz條件，可用數(shù)學(xué)歸納法證明所有的關(guān)于滿足Lipschitz條件，從而增量函數(shù)關(guān)于滿足Lipschitz條件，所以Runge-Kutta方法收斂。而且對于這類滿足Lipschitz條件的單步法而言，相容條件是收斂條件的充分必要條件。</p><p><b>　　III、穩(wěn)定性</b></p><p>　　下

119、面討論穩(wěn)定性問題，更準(zhǔn)確地說就是數(shù)值解關(guān)于初值的連續(xù)依賴性問題。</p><p>　　定義4 如果一個單步法滿足：存在正常數(shù)和，使得用任意步長以及任意兩個初值和計算得到的兩組數(shù)值和成立</p><p><b>　　，</b></p><p>　　則稱該方法是穩(wěn)定的，其中。</p><p>　　定理2 若增量函數(shù)關(guān)于滿

120、足Lipschitz條件，則公式（2）是穩(wěn)定的。</p><p>　　由上可知，歐拉法，改進(jìn)歐拉法和Runge-Kutta方法都是穩(wěn)定的。但是因為上面的穩(wěn)定性定義中沒有考慮計算過程中的舍入誤差，使得最后得到的結(jié)果可能會很差，所以又有了任何一步產(chǎn)生的誤差在以后過程中足部減少的絕對穩(wěn)定性。它不僅與函數(shù)有關(guān)，還和步長有關(guān)。</p><p><b>　　對于模型問題：</b>

121、</p><p><b>　　， </b></p><p><b>　?。?0）</b></p><p>　　其中是常數(shù)（可以是復(fù)數(shù)）。</p><p>　　定義5 對初值問題（10），記</p><p><b>　　，</b></p>

122、;<p>　　若一個數(shù)值方法在某一步由計算產(chǎn)生的誤差在以后各步中逐步減少，則稱該方法關(guān)于是絕對穩(wěn)定的。復(fù)平面上所有這樣的組成的區(qū)域稱為該方法的絕對穩(wěn)定區(qū)域。</p><p>　　由此可知，當(dāng)滿足時，n級n階的Runge-Kutta方法(n=1,2,3,4)是絕對穩(wěn)定的。</p><p>　　如果一個方法的絕對穩(wěn)定區(qū)域比較大，則可說它的絕對穩(wěn)定性比較好，且隱式方法的絕對穩(wěn)定性比

123、顯式方法好；對于Runge-Kutta方法，階愈高絕對穩(wěn)定性愈好。</p><p><b>　　三、總結(jié)部分 </b></p><p>　　常微分方程在科學(xué)和工程中常用于建立數(shù)學(xué)模型，通常它們沒有解析解，而需要數(shù)值方法來近似求解。在科學(xué)的計算機(jī)化進(jìn)程中，常微分方程也在不斷地精進(jìn)中，更在不同的領(lǐng)域得到了前所未有的發(fā)展和應(yīng)用。求解常微分方程的初值問題的數(shù)值方法有單步法和

124、多步法。其中單步法中的Runge-Kutta方法是目前工程上求常微分方程數(shù)值解的基本方法。但是這并不表示Runge-Kutta方法是最好的，它也存在著缺點，例如：它計算函數(shù)值的次數(shù)較多。尤其是在級數(shù)高于4級時，這一缺點特別明顯：計算函數(shù)值的工作量增加較快，但是精度即階數(shù)提高較慢。還有各階Runge-Kutta方法的絕對穩(wěn)定區(qū)域不夠大，對于壞條件的初值問題不盡適用等。但是它作為一類便于應(yīng)用的單步方法，其實用價值還是很高的，特別是在現(xiàn)如今計