

版權說明:本文檔由用戶提供并上傳,收益歸屬內容提供方,若內容存在侵權,請進行舉報或認領
文檔簡介
1、<p> 本科畢業(yè)論文(設計)</p><p><b> ( 20 屆)</b></p><p> 常微分方程初值問題的Runge-Kutta解法</p><p> 所在學院 </p><p> 專業(yè)班級 信息與計算科學
2、 </p><p> 學生姓名 學號 </p><p> 指導教師 職稱 </p><p> 完成日期 年 月 </p><p> 摘要:常微分方程是指只有一個自變量的微分方程。常微分方程在很多學
3、科領域內有著重要的作用,如自動控制、各種電子學裝置的設計等。本文首先介紹了常微分方程以及初值問題的除Runge-Kutta解法外的數值解一些理論;接著,我們對常微分方程初值問題提出Runge-Kutta解法的構造公式;然后,我們再對這個公式進行誤差分析,包括相容性、穩(wěn)定性和收斂性分析;最后,通過一些數值實例,我們驗證了Runge-Kutta法的優(yōu)缺點。</p><p> 關鍵詞:常微分方程;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 常微分方程初值問題的一些數值解法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解法的構造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 數值算例14</p><p
10、> 5 結 論16</p><p><b> 致 謝17</b></p><p><b> 參考文獻18</b></p><p><b> 附錄19</b></p><p><b> 1 緒 論</b></p>
11、<p> 1.1 問題的背景</p><p> 常微分方程在微積分概念出現(xiàn)后即已出現(xiàn)。從萊布尼茲專門研究用變量變換解決一階微分方程的求解問題的“求通解”時代,到1841年劉維爾的里卡蒂方程不存在一般初等解和柯西的初值問題的提出,常微分方程轉向“求定解”時代。再到19世紀末天體力學的太陽系穩(wěn)定性問題眼界需要常微分方程解的大范圍性態(tài),從而發(fā)展到了“求所有解問題”。最后到了20世紀六七十年代,常微分方
12、程因為計算機的發(fā)展進入“求特殊解”階段,發(fā)現(xiàn)了具有新性質的特殊的解和方程,如混沌(解)、奇異吸引子及孤立子等等[1]。</p><p> 常微分方程是指只有一個自變量的微分方程。而且常微分方程在很多學科領域內有著重要的作用,如自動控制、各種電子學裝置的設計、彈道的計算、飛機和導彈飛行的穩(wěn)定性的研究、化學反應過程穩(wěn)定性的研究等等。常微分方程初值問題的數值求解是求近似解的一種經典方法。</p><
13、;p> 1.2 常微分方程初值問題簡介[2]-[9]</p><p> 大多數學科領域內的問題需要求解一個初值問題,即求解滿足給定初值條件的微分方程。我們更多的是使用逼近原問題的解的方法來逼近原方程的解。因逼近方法給出更精確的結果和實際的誤差信息,所以更經常被使用。現(xiàn)代科學、技術、工程中的大量數學模型都可以用常微分方程的初值問題來描述,所以研究常微分方程初值問題的一些性質、解法很有必要。一些典型的常微
14、分初值問題的數值求解問題的方法有:單步法和線性多步法。</p><p> 在求解區(qū)間上取定節(jié)點</p><p> 令,稱為積分網格的步長。用表示初值問題精確解在節(jié)點上函數值的近似值。對給定的數值積分方法,中的各個市按某一個遞推算法確定的。一個遞推算法,如果在用它計算時只用到已經求出的諸值中的,而無須使用其余值中的任何一個,則稱此算法為單步方法。相反的則是多步法。</p>
15、<p> 對于一階的常微分方程初值問題</p><p><b> , </b></p><p><b> (1)</b></p><p> 的數值方法,其中是和的已知函數,為給定的初值。為使上面的解存在、唯一且連續(xù)依賴初值,則必須關于滿足Lipschitz條件:即存在常數L,使得對于任意,均成立不
16、等式</p><p> 單步顯式公式的一般形式為</p><p><b> ,</b></p><p><b> ?。?)</b></p><p> 稱為增量函數。一般來說微分方程問題(1)的精確解不滿足式(2),即一般地有</p><p><b> 定義1
17、 稱</b></p><p> 為單步顯式公式(2)在處的截斷誤差。</p><p> 一個求解公式的局部截斷誤差刻畫了其逼近微分方程的準確程度。根據上述定義可以直接求得各式的局部截斷誤差。</p><p> 定義2 如果一個求解公式的局部截斷誤差為,則稱該求解公式是階的,或稱具有階精度。</p><p> 在多步
18、法中也有幾個基本的理論問題的定義,定理。</p><p> 定義3 對于初值問題(1)的多步法</p><p> , (3)</p><p> 說是相容的,如果它至少是一階的。</p><p> 定理1 為使步法(3)相容,必需且只需</p><p>
19、 2 常微分方程初值問題的一些數值解法</p><p> 下面將先簡單的介紹一下除Runge-Kutta法外的一些數值方法。</p><p> 2.1 單步法中的歐拉解法[4]-[9]</p><p> 單步法中最簡單的數值解法是歐拉法。歐拉法的基本思想是由推導出下一個節(jié)點,其推導公式為</p><p><b> ,.&
20、lt;/b></p><p> 其局部截斷誤差為。歐拉法的幾何意義是用一條過的折線來近似替代過的積分曲線。</p><p> 雖然這個方法的計算很簡單,但精確度較低,所以有了精確度較高的梯形方法。推導公式為,,它的局部截斷誤差為。此為隱型格式,迭代求解所需計算量大。</p><p> 結合上述兩種方法的優(yōu)點,從而得到改進歐拉方法。其計算公式為, 。
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解法的構造[5]-[15]</
23、p><p> 用Taylor級數法的公式,可以構造高階單步法,但增量函數是用的各階導數表示的,通常不易計算。由于函數在一點的導數值可以用該點附近若干點的函數值近似表示,因此可以把Taylor級數法中的增量函數改為在一些點上函數值的組合,然后用Taylor展開確定待定的系數,使方法達到一定的階。這就是Runge-Kutta法的基本思想?;蛘呤窃诿總€區(qū)間多取幾點,將它們的斜率加權平均作為平均斜率,就有可能構造出精度更高
24、的計算公式。通過構造高階單步法來提高精度,而較高的精度意味著計算結果更加精確,誤差會隨著的減小而迅速減小。</p><p><b> 現(xiàn)在假設</b></p><p><b> ?。?)</b></p><p><b> 其中</b></p><p><b>
25、,</b></p><p> , (5)</p><p> 此類遞推算法為顯示N級Runge-Kutta方法。上述算法公式中的系數,和希望如此確定,使得(4)式的Taylor展開式中所有冪次不超過的那些項的系數與Taylor展開式 (6) (其中,)中的相應項的系數相等。</p&g
26、t;<p> 當N=1時,則可知=1,此時式(4)為歐拉公式。</p><p> 當N=2時,確定系數,和使得與兩者的Taylor展開中含項的系數,,對盡可能大的相等。</p><p> 利用二元的Taylor展開公式</p><p><b> 可將(4)式表示為</b></p><p><b
27、> (7)</b></p><p> 再將Taylor展式(4)改成</p><p><b> ?。?)</b></p><p> 比較(7)和(8)式中項的系數,得到</p><p> , , (9)</p><
28、p> 因為有4個待定系數,可是就只有3個方程數目。因此具有無窮多個數組滿足條件(9),即存在無窮多個二級Runge-Kutta方法,它們的截斷誤差。</p><p><b> 若取,便得到</b></p><p><b> 或</b></p><p><b> 這是改進歐拉公式。</b>
29、</p><p> 因為二級的RK方法最多是二階的,若想得到更高階的方法,必須增加級數,即取較大的N。當N較大時,計算也就會更加的復雜。</p><p> 當N=3時,和N=2的推導過程相相似,當時,得到</p><p><b> 或當時,有</b></p><p> 同理,N=4時,構造經典的四階RK方法,其計
30、算公式為</p><p><b> ?。?0)</b></p><p><b> 其局部截斷誤差為。</b></p><p> 但是,以上的方法還是不夠的。在涉及初值問題數值解的實際問題中,總要先設定一個精度允許值,數值解與真解的偏離不能超過這個允許值。因此,步長還受到這個允許值的控制,常常會發(fā)生這樣的情形:僅在求解的
31、某一部分需要小的步長,而在其余部分則用較大的步長就可以了。</p><p> 所以,各種用于算法中自動地調整步長的方法(稱為自適應方法)就得到了發(fā)展。就如在四階的Runge-Kutta方法中,為了從到,我們可以用步長為做一步,也可以用步長為做兩步,如果沒有截斷誤差,它們求出的的值應是相同的。但事實上它們不會相同,所以其差異可看成是對局部截斷誤差的一個估計:如果這個差異不超過預先給定的允許值,那么步長就是滿意的;
32、如果它超過允許值,就將減半;如果它遠小于允許值,就將步長加倍。</p><p> 雖然這樣的過程很容易計算,但要浪費大量的時間,因而不宜提倡。</p><p> 所以給出一個Runge-Kutta方法的漂亮的自適應格式:</p><p><b> ?。?1)</b></p><p> 它要求計算五個函數值,但只是一
33、個四階方法,因此單獨使用并無好處,但它可以通過增加一個函數值(而不需要改變已有的五個函數值)來得到一個更高一階的方法:</p><p> 同(11)式 (12)</p><p> 從四階與五階過程中得到的的值之間的差,是對四階方法中局部截斷誤差的一個估計。換句話說,通過求六個函數值和簡單的四則運算(系數可以預先計算
34、好存放起來),我們得到了一個四階的近似值連同一個誤差估計!計算機只要比較這個近似值和允許誤差界,就可以自動地決定在下一步計算中,到底是保持還是縮放當前正在使用的步長。</p><p> 由上可知,由于RK方法是通過對函數的Taylor展開來構造的,因此高階RK方法要求微分方程的解具有高階導數。如果解的光滑性較差,即使用高階方法也不可能得到高精度的數值解,此時 應當用低階方法?;蛘哒f從這些上述公式的討論中可知,級
35、數N表示計算函數值的次數,當時,N級的RK方法的階數,即計算函數值的次數越多,RK方法的精度即階數越高。但是當時,情況就大不一樣了,J.C.Butcher指出:當時,;當時,;而且,。由此可見四級以上的RK方法,計算函數值的工作量增加較快,而階數即精度提高較慢,故而,在實際問題中比較常用的是四級RK方法。</p><p> 3.2 誤差分析[9]-[15]</p><p> 3.2.
36、1 相容性</p><p> 定義4 當,時,稱數值方法是相容的。</p><p><b> 由于</b></p><p><b> 于是</b></p><p> 故有結論:R-K法相容。即p階的Runge-Kutta方法是相容的。</p><p> 3.2.
37、2 收斂性</p><p> 由于截斷誤差的存在,計數值只是精確解的一個近似。故我們可以相當于有這樣的問題:差是不是隨著步長的減小而減小,這就是收斂性問題。</p><p> 定義5 當時,,,稱數值方法是收斂的。</p><p><b> 令 ,則</b></p><p><b> -) <
38、;/b></p><p><b> 于是得到結論:</b></p><p> 當時,,且有誤差估計</p><p><b> 。</b></p><p> 即也有由于,當時,,當步長減小時,從計算到點的步數就增大,從而推出下面的定理。</p><p> 定理2
39、 若公式(2)是階方法,且增量函數關于滿足Lipschitz條件,則該方法收斂。</p><p> 對于Runge-Kutta方法,當函數關于滿足Lipschitz條件,可用數學歸納法證明所有的關于滿足Lipschitz條件,從而增量函數關于滿足Lipschitz條件,所以Runge-Kutta方法收斂。而且對于這類滿足Lipschitz條件的單步法而言,相容條件是收斂條件的充分必要條件。</p>
40、<p> 3.2.3 穩(wěn)定性</p><p> 下面討論穩(wěn)定性問題,更準確地說就是數值解關于初值的連續(xù)依賴性問題。</p><p> 定義6 如果一個單步法滿足:存在正常數和,使得用任意步長以及任意兩個初值和計算得到的兩組數值和成立</p><p><b> ,</b></p><p> 則稱
41、該方法是穩(wěn)定的,其中。</p><p> 定理3 若增量函數關于滿足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、以上是在的條件下才滿足的。故可知,歐拉法,改進歐拉法和Runge-Kutta方法都是穩(wěn)定的。然而,第一,我們實際上是取有限的固定步長進行計算的,它并不能隨意縮小至與零充分接近。第二,即使可以任意縮小,我們也不會那么做,這將導致計算工作亮的急劇增加。第三,更重要的是,由于舍入誤差在計算過程中的產生、傳播和放大,小的步長可能導致計算結果嚴重失真甚至面目全非。因此在現(xiàn)實中我們關心的恰恰是對一個方法來講,它允許的最大步長是多少。顯然,這個最大步
44、長與所給方程的性質有關。所以就又有了使任何一步產生的誤差在以后過程中足夠減少的絕對穩(wěn)定性。</p><p> 定義7 對于Runge-Kutta方法的</p><p><b> 當即充分小時,</b></p><p> Runge-Kutta方法穩(wěn)定滿足根條件</p><p><b> 當不趨近0時,
45、</b></p><p><b> -) </b></p><p><b> 一般地,記 ,,則</b></p><p> 其中是的次齊次多項式,且其系數之和為。于是</p><p><b> 即</b></p><p> Run
46、ge-Kutta方法絕對穩(wěn)定滿足嚴格根條件。即求,使</p><p><b> 滿足:</b></p><p> 這時,所求的的變動范圍就是絕對穩(wěn)定區(qū)域。</p><p> 下表3-1給出了高階單步方法在實數域中的絕對穩(wěn)定區(qū)間:</p><p> 表3-1絕對穩(wěn)定區(qū)間</p><p>
47、由上表可以發(fā)現(xiàn),當增加時它的絕對穩(wěn)定區(qū)域還有所增大[15]。</p><p> 故由此可知,當滿足時,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 數值算例</b></p><p> 因為現(xiàn)在無法比較出4階的Runge-Kutta方法與其他階數求解方法的優(yōu)劣,所以這里引進比較低階Runge-Kutta方法的一個度量描述如下:</p><p> 4階Runge-Kutta方法每步需要4次求值,所以如果它有更優(yōu)的話,它應能比1/4步長的歐拉法給出更精確地答案。類似地,如果4階Runge-Kutta方
49、法比二階Runge-Kutta方法優(yōu)越的話,則步長為的這種方法應該比步長為的二階方法給出更精確地結果,因為4階方法每步需要二倍次數的求值。</p><p> 算例4-1 對于問題</p><p><b> ,,,</b></p><p> 用四階的Runge-Kutta方法解上式問題。</p><p> 解:對
50、于上式我們采用不同階數的Runge-Kutta方法,即,,,時,令步長,從而得出各種不同階數的Runge-Kutta方法的相對誤差,如下表4-1所示。</p><p><b> 表4-1 </b></p><p> 由以上的表可以觀察到4階Runge-Kutta方法的誤差是最小的。即比其它的階數Runge-Kutta方法有更高的精確度。</p>&l
51、t;p> 算例4-2 對于問題</p><p><b> ,,</b></p><p> 分別取步長為(i=1,2),求解。</p><p> 解:當時,不屬于絕對穩(wěn)定區(qū)域,而步長時,相應的屬于絕對穩(wěn)定區(qū)域</p><p> 計算結果如下表4-2所示。</p><p><b
52、> 表4-2</b></p><p> 從上表中,我們可以觀察到當時,其相對誤差是呈現(xiàn)出收斂與常數的趨勢。但當時,其相對誤差就呈現(xiàn)發(fā)散的趨勢,顯然這是由于此時的不屬于絕對穩(wěn)定區(qū)域所引起的,故我們可以知道絕對穩(wěn)定性的重要。</p><p><b> 5 結 論</b></p><p> 本文先介紹了常微分方程的初值問題
53、的一些數值解法,在詳細介紹了單步法中的Runge-Kutta解法。常微分方程是指只有一個自變量的微分方程,而且常微分方程在很多學科領域內有著重要的作用。常微分方程初值問題的數值求解是求近似解的一種經典方法,Runge-Kutta方法是目前工程上求常微分方程數值解的基本方法。它的基本思想是在每個區(qū)間多取幾點,將它們的斜率加權平均作為平均斜率,就有可能構造出精度更高的計算公式。</p><p> 通過對Runge-
54、Kutta方法的構造,通過對Runge-Kutta方法的研究,得出了Runge-Kutta法的相容性、收斂性和穩(wěn)定性分析,通過程序(MATLAB語言)編程,給出數值例子,尤其是4階的Runge-Kutta公式,驗證了理論分析。常微分方程初值問題的數值解在數學意義上的存在性可以在一般的條件下得到證明,而從實際應用上來講,人們需要的往往并不是解在數學上的存在性,而是在我們所關心的某個定義范圍內,對應于某些特定的自變量的取值,尋求數值解,而這
55、也是為什么要對其數值計算法中Runge-Kutta法的研究。</p><p><b> 參考文獻</b></p><p> [1] 王高雄,朱思銘,王壽松.常微分方程[M].北京:高等教育出版社,2006.</p><p> [2] [美]Richard L.Burden,J.Douglas Faires. 數值分析[[M].北京:高等教
56、育出版社,2005.</p><p> [3] 中國大百科全書·數學.北京:中國大百科全書出版社,1988.</p><p> [4] 余徳浩,湯華中.微分方程數值解法[M].北京:科學出版社,2003.</p><p> [5] 胡健偉,湯懷民.微分方程數值解法[M].北京:科學出版社,1999.</p><p> [6]
57、 李榮華,馮果忱.微分方程數值解法[M].北京:高等教育出版社,1997.</p><p> [7] Arieh Iserles.微分方程數值分析基礎教程[M].北京:清華大學出版社,2005.</p><p> [8] 李忠杰.常微分方程初值問題RK法和多步法[J].黑龍江科技信息,科教文化,(1):199.</p><p> [9] 李瑞遐,何志慶.微分方
58、程數值方法[M].上海:華東理工大學出版社,2005.</p><p> [10] 孫維夫,姜云義.常微分方程初值問題RK法和多步法[J].煙臺職業(yè)學院學報,2010,(2):84-88.</p><p> [11] 劉利斌,王偉.四階Runge-Kutta算法的優(yōu)化分析[J].成都大學學報(自然科學版),2007,(1):19-21.</p><p> [1
59、2] 林群.微分方程數值解法基礎教程[M].北京:科學出版社,2003.</p><p> [13] 馬明書.龍格-庫塔法的級數和階數[J].河南師范大學河南機專學報,1994,(2):5-7.</p><p> [14] 吳怡,張向輝,熊超西.三階Runge-Kutta最優(yōu)算法[J].河北秦皇島職業(yè)技術學院基礎部,2005,(18):57-61.</p><p&g
60、t; [15] 李立康,於崇華.微分方程數值方法[M].上海:復旦大學出版社,2003.</p><p><b> 附錄</b></p><p> 附錄1:例題4.1的編程</p><p><b> 函數的表達式:</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表示函數在x0處的值,可以是向量</p><p> %PointNum表示自變量在[x0,xt]上取的點數,包括端點.</p><p&g
66、t; if nargin<5 | PointNum<=0 %如果函數僅輸入4個參數值,則PointNum默認值為101</p><p> PointNum=101;</p><p><b> end</b></p><p> if nargin<4 %y0默認值為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;%自變量數組</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)%畫出方程解的函數圖</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為離散點數,varargin為可選輸入項可傳適當參數給函數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為離散點數,varargin為可選輸入項可傳適當參數給函數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為離散點數,varargin為可選輸入項可傳適當參數給函數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、> 文獻綜述</b></p><p> 常微分方程初值問題的Runge-Kutta解法 </p><p><b> 一、前言部分</b></p><p> 常微分方程在很多學科領域內有著重要的作用,如自動控制、各種電子學裝置的設計、彈道的計算、飛機和導彈飛行的穩(wěn)定性的研究、化學反應
94、過程穩(wěn)定性的研究等等,這些問題都可以化為求微分方程的解,或者化為研究解的性質的問題。這些問題都包含某個變量關于另一個變量的變化。大多數這樣的問題需要求解一個初值問題,即求解滿足給定初值條件的微分方程。我們更多的是使用逼近原問題的解的方法來逼近原方程的解。因逼近方法給出更精確的結果和實際的誤差信息,所以更經常被使用。一些典型的常微分初值問題的數值求解問題的方法有:單步法和線性多步法。</p><p> 在求解區(qū)間
95、上取定節(jié)點</p><p> 令,稱為積分網格的步長。用表示初值問題精確解在節(jié)點上函數值的近似值。對給定的數值積分方法,中的各個市按某一個遞推算法確定的。一個遞推算法,如果在用它計算時只用到已經求出的諸值中的,而無須使用其余值中的任何一個,則稱此算法為單步方法。相反的則是多步法。單步法主要有歐拉法和Runge-Kutta法,多步法主要有Adams法和Milne法等。</p><p>
96、本文綜述常微分初值問題初值問題的數值解法及其誤差估計(相容性、穩(wěn)定性和收斂性分析),重點介紹了Runge-Kutta法。</p><p><b> 二、主題部分</b></p><p> 2.1 常微分方程的初值問題概述</p><p> 常微分方程在微積分概念出現(xiàn)后即已出現(xiàn)。從萊布尼茲專門研究用變量變換解決一階微分方程的求解問題的“求通
97、解”時代,到1841年劉維爾的里卡蒂方程不存在一般初等解和柯西的初值問題的提出,常微分方程轉向“求定解”時代。再到19世紀末天體力學的太陽系穩(wěn)定性問題眼界需要常微分方程解的大范圍性態(tài),從而發(fā)展到了“求所有解問題”。最后到了20世紀六七十年代,常微分方程因為計算機的發(fā)展進入“求特殊解”階段,發(fā)現(xiàn)了具有新性質的特殊的解和方程,如混沌(解)、奇異吸引子及孤立子等等。</p><p> 在數學分析(微積分)中研究了變量
98、的各種函數及函數的微分與積分。如函數未知,但知道變量與函數的代數關系式,便組成代數方程,通過求解代數方程解出未知函數。同樣,如果知道自變量、未知函數及函數的導數(或微分)組成的關系式,得到的便是微分方程,通過求解微分方程求出未知函數。自變量只有一個的微分方程稱為常微分方程。常微分方程是數學分析或基礎數學的一個組成部分,在整個數學大廈中占據著重要位置。</p><p> 在常微分方程中,我們知道只有少數簡單的常微
99、分方程才能用初等積分法求出它們的精確解析解,多數情形只能用近似解法求其近似解。近似解法可以分成兩類:解析方法和數值方法。其中前者尋求解的近似表達式,后者則是計算微分方程解在求解區(qū)域中一些離散點上的近似值。</p><p> 對于一階的常微分方程初值問題</p><p><b> , </b></p><p><b> ?。?/p>
100、1)</b></p><p> 的數值方法,其中是和的已知函數,為給定的初值。為使上面的解存在、唯一且連續(xù)依賴初值,則必須關于滿足Lipschitz條件:即存在常數L,使得對于任意,均成立不等式</p><p> 單步顯式公式的一般形式為</p><p><b> ,</b></p><p><b
101、> ?。?)</b></p><p> 稱為增量函數。一般來說微分方程問題(1)的精確解不滿足式(2),即一般地有</p><p><b> 定義1 稱</b></p><p> 為單步顯式公式(2)在處的截斷誤差。</p><p> 一個求解公式的局部截斷誤差刻畫了其逼近微分方程的準確程度
102、。根據上述定義可以直接求的各式的局部截斷誤差。例如歐拉公式的局部截斷誤差為</p><p> 定義2 如果一個求解公式的局部截斷誤差為,則稱該求解公式是階的,或稱具有階精度。</p><p> 2.2 常微分方程的初值問題的數值方法</p><p> 下面將先簡單的介紹一下除Runge-Kutta法外的一些數值方法。</p><p>
103、; 2.2.1 單步法中的歐拉法</p><p> 單步法中最簡單的數值解法是歐拉法。歐拉法的基本思想是由推導出下一個節(jié)點,其推導公式為</p><p><b> ,.</b></p><p> 其局部截斷誤差為。歐拉法的幾何意義是用一條過的折線來近似替代過的積分曲線。</p><p> 雖然這個方法的計算很
104、簡單,但精確度較低,所以有了精確度較高的梯形方法。推導公式為,,它的局部截斷誤差為。此為隱型格式,迭代求解所需計算量大。</p><p> 結合上述兩種方法的優(yōu)點,從而得到改進歐拉方法。其計算公式為, 。</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法的構造</p><p> 用Taylor級數法的公式,可以構造高階單步法,但增量函數是用的各階導數表示的,通常不易計算。由于函數在一點的導數值可以用該點附近若干點的函數值近似表示,因此可以把Taylor級數法中的增量
107、函數改為在一些點上函數值的組合,然后用Taylor展開確定待定的系數,使方法達到一定的階。這就是Runge-Kutta法的基本思想?;蛘呤窃诿總€區(qū)間多取幾點,將它們的斜率加權平均作為平均斜率,就有可能構造出精度更高的計算公式。通過構造高階單步法來提高精度,而較高的精度意味著計算結果更加精確,誤差會隨著的減小而迅速減小。</p><p><b> 現(xiàn)在假設</b></p>&l
108、t;p><b> ?。?)</b></p><p><b> 其中</b></p><p><b> ,</b></p><p> , (4)</p><p> 此類遞推算法為顯示n級Runge-Kutta方法。上述算法公式中的系數,和希望如此確定,使得
109、(3)式的Taylor展開式中所有冪次不超過的那些項的系數與Taylor展開式 (5)(其中,)中的相應項的系數相等。</p><p> 當n=1時,則可知=1,此時式(3)為歐拉公式。</p><p> 當n=2時,確定系數,和使得與兩者的Taylor展開中含項的系數,,對盡可能大的相等。</p><p> 利用二元的Taylor
110、展開公式</p><p><b> 可將(3)式表示為</b></p><p><b> ?。?)</b></p><p> 再將Taylor展式(3)改成</p><p><b> ?。?)</b></p><p> 比較(6)和(7)式中項的系
111、數,得到</p><p> , , (8)</p><p> 因為有4個待定系數,可是就只有3個方程數目。因此具有無窮多個數組滿足條件(8),即存在無窮多個二級Runge-Kutta方法,它們的截斷誤差。</p><p><b> 若取,便得到</b></p><p&
112、gt;<b> 或</b></p><p><b> 這是改進歐拉公式。</b></p><p> 因為二級的RK方法最多是二階的,若想得到更高階的方法,必須增加級數,即取較大的n。當n較大時,計算也就會更加的復雜。</p><p> 當n=3時,和n=2的推導過程相相似,當時,得到</p><
113、p><b> 或當時,有</b></p><p> 同理,n=4時,構造經典的四階RK方法,其計算公式為</p><p><b> ?。?)</b></p><p><b> 其局部截斷誤差為。</b></p><p> 由于RK方法是通過對函數的Taylor展開
114、來構造的,因此高階RK方法要求微分方程的解具有高階導數。如果解的光滑性較差,即使用高階方法也不可能得到高精度的數值解,此時 應當用低階方法?;蛘哒f從這些上述公式的討論中可知,級數n表示計算函數值的次數,當時,n級的RK方法的階數,即計算函數值的次數越多,RK方法的精度即階數越高。但是當時,情況就大不一樣了,J.C.Butcher指出:當時,;當時,;而且,。由此可見四級以上的RK方法,計算函數值的工作量增加較快,而階數即精度提高較慢,故
115、而,在實際問題中比較常用的是四級RK方法。</p><p> 2.3.2 Runge-Kutta法的誤差分析</p><p><b> I、相容性</b></p><p> 當,時,稱數值方法是相容的。</p><p><b> 由于</b></p><p><
116、;b> 于是</b></p><p> 故有結論:R-K法相容。即p階的Runge-Kutta方法是相容的。</p><p><b> II、收斂性</b></p><p> 由于截斷誤差的存在,計數值只是精確解的一個近似。故我們可以相當這樣的問題:差是不是隨著步長的減小而減小。這就是收斂性問題。</p>
117、<p> 定義3 如果一個單步法滿足:對中任意一個固定的點成立</p><p> 則稱該方法是收斂的。</p><p> 由于,當時,,也即當步長減小時,從計算到點的步數就增大,故有下面的定理。</p><p> 定理1 若公式(2)是階方法,且增量函數關于滿足Lipschitz條件,則該方法收斂。</p><p>
118、 對于Runge-Kutta方法,當函數關于滿足Lipschitz條件,可用數學歸納法證明所有的關于滿足Lipschitz條件,從而增量函數關于滿足Lipschitz條件,所以Runge-Kutta方法收斂。而且對于這類滿足Lipschitz條件的單步法而言,相容條件是收斂條件的充分必要條件。</p><p><b> III、穩(wěn)定性</b></p><p> 下
119、面討論穩(wěn)定性問題,更準確地說就是數值解關于初值的連續(xù)依賴性問題。</p><p> 定義4 如果一個單步法滿足:存在正常數和,使得用任意步長以及任意兩個初值和計算得到的兩組數值和成立</p><p><b> ,</b></p><p> 則稱該方法是穩(wěn)定的,其中。</p><p> 定理2 若增量函數關于滿
120、足Lipschitz條件,則公式(2)是穩(wěn)定的。</p><p> 由上可知,歐拉法,改進歐拉法和Runge-Kutta方法都是穩(wěn)定的。但是因為上面的穩(wěn)定性定義中沒有考慮計算過程中的舍入誤差,使得最后得到的結果可能會很差,所以又有了任何一步產生的誤差在以后過程中足部減少的絕對穩(wěn)定性。它不僅與函數有關,還和步長有關。</p><p><b> 對于模型問題:</b>
121、</p><p><b> , </b></p><p><b> ?。?0)</b></p><p> 其中是常數(可以是復數)。</p><p> 定義5 對初值問題(10),記</p><p><b> ,</b></p>
122、;<p> 若一個數值方法在某一步由計算產生的誤差在以后各步中逐步減少,則稱該方法關于是絕對穩(wěn)定的。復平面上所有這樣的組成的區(qū)域稱為該方法的絕對穩(wěn)定區(qū)域。</p><p> 由此可知,當滿足時,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> 三、總結部分 </b></p><p> 常微分方程在科學和工程中常用于建立數學模型,通常它們沒有解析解,而需要數值方法來近似求解。在科學的計算機化進程中,常微分方程也在不斷地精進中,更在不同的領域得到了前所未有的發(fā)展和應用。求解常微分方程的初值問題的數值方法有單步法和
124、多步法。其中單步法中的Runge-Kutta方法是目前工程上求常微分方程數值解的基本方法。但是這并不表示Runge-Kutta方法是最好的,它也存在著缺點,例如:它計算函數值的次數較多。尤其是在級數高于4級時,這一缺點特別明顯:計算函數值的工作量增加較快,但是精度即階數提高較慢。還有各階Runge-Kutta方法的絕對穩(wěn)定區(qū)域不夠大,對于壞條件的初值問題不盡適用等。但是它作為一類便于應用的單步方法,其實用價值還是很高的,特別是在現(xiàn)如今計
溫馨提示
- 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
- 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權益歸上傳用戶所有。
- 3. 本站RAR壓縮包中若帶圖紙,網頁內容里面會有圖紙預覽,若沒有圖紙預覽就沒有圖紙。
- 4. 未經權益所有人同意不得將文件中的內容挪作商業(yè)或盈利用途。
- 5. 眾賞文庫僅提供信息存儲空間,僅對用戶上傳內容的表現(xiàn)方式做保護處理,對用戶上傳分享的文檔內容本身不做任何修改或編輯,并不能對任何下載內容負責。
- 6. 下載文件中如有侵權或不適當內容,請與我們聯(lián)系,我們立即糾正。
- 7. 本站不保證下載資源的準確性、安全性和完整性, 同時也不承擔用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。
最新文檔
- 常微分方程初值問題的runge-kutta解法[文獻綜述]
- 常微分方程初值問題的runge-kutta解法[開題報告]
- 常微分方程初值問題的runge-kutta解法[畢業(yè)論文]
- 常微分方程初值問題的預估-校正解法【信息科學與技術專業(yè)】【畢業(yè)設計+文獻綜述+開題報告】
- 常微分方程初值問題的預估-校正解法[文獻綜述]
- 淺談常微分方程的數值解法及其應用【信息科學與技術專業(yè)】【畢業(yè)設計+文獻綜述+開題報告】
- 常微分方程初值問題的預估-校正解法[開題報告]
- 常微分方程初值問題的預估-校正解法[畢業(yè)論文]
- 解常微分方程的三步Runge-Kutta方法.pdf
- 矩陣方程的數值解法【信息科學與技術專業(yè)】【畢業(yè)設計+文獻綜述+開題報告】
- 求解振蕩常微分方程的辛指數擬合Runge-Kutta(-Nystrom)型方法.pdf
- 求解分數階微分方程的Runge-Kutta方法.pdf
- 剛性延遲積分微分方程的Runge-Kutta離散.pdf
- 幾類常微分方程典型的解法[畢業(yè)論文+開題報告+文獻綜述]
- 延遲積分微分方程的變步長Runge-Kutta方法.pdf
- 橢圓型偏微分方程的求解及其應用【信息科學與技術專業(yè)】【畢業(yè)設計+文獻綜述+開題報告】
- 雙曲型偏微分方程的求解及其應用【信息科學與技術專業(yè)】【畢業(yè)設計+文獻綜述+開題報告】
- 幾類求解分數階微分方程的Runge-Kutta方法.pdf
- 模糊微分方程的初值問題.pdf
- 延遲積分微分方程多步Runge-Kutta法及并行實現(xiàn).pdf
評論
0/150
提交評論