MATLAB第十三講常微分方程初值問題數(shù)值解法

1、上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 常微分方程初值問常微分方程初值問題數(shù)值解法題數(shù)值解法上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)9.1 引引 言言 科學(xué)技術(shù)中很多問題都可用微分方程的科學(xué)技術(shù)中很多問題都可用微分方程的定解問定解問題題來(lái)描述，主要有來(lái)描述，主要有初值問題初值問題與與邊值問題邊值問題兩大類，兩大類，本章本章只考慮初值問題只考慮初值問題. 常微分方程初值問題中最簡(jiǎn)單的例常微分方程初值問題中最簡(jiǎn)單的例子是人口模型，設(shè)某特定區(qū)域在子是人口模型，設(shè)某特定區(qū)域在t0時(shí)刻人口為時(shí)刻人口為y(t0)=y0已知的，該區(qū)域的人口自然增長(zhǎng)率為已知的，該區(qū)域的人口自然增長(zhǎng)率為 ，人

2、口增長(zhǎng)與，人口增長(zhǎng)與人口總數(shù)成正比，所以人口總數(shù)成正比，所以t時(shí)刻的人口總數(shù)時(shí)刻的人口總數(shù)y(t)滿足以下滿足以下微分方程微分方程 .)(),(00ytytyy 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)很多物理系統(tǒng)與時(shí)間有關(guān)，從衛(wèi)星運(yùn)行軌道到單很多物理系統(tǒng)與時(shí)間有關(guān)，從衛(wèi)星運(yùn)行軌道到單擺運(yùn)動(dòng)，從化學(xué)反應(yīng)到物種競(jìng)爭(zhēng)都是隨時(shí)間的延續(xù)而擺運(yùn)動(dòng)，從化學(xué)反應(yīng)到物種競(jìng)爭(zhēng)都是隨時(shí)間的延續(xù)而不斷變化的不斷變化的. 解常微分方程是描述連續(xù)變化的數(shù)學(xué)語(yǔ)解常微分方程是描述連續(xù)變化的數(shù)學(xué)語(yǔ)言，微分方程的求解就是確定滿足給定方程的可微函言，微分方程的求解就是確定滿足給定方程的可微函數(shù)數(shù)y(t)，研究它的數(shù)值方法

3、是本章的主要目的，研究它的數(shù)值方法是本章的主要目的. 考慮考慮一一階常微分方程的初值問題階常微分方程的初值問題 )2 . 1(.)()1 . 1(,),(000yxybxxyxfy則稱則稱f關(guān)于關(guān)于y滿足滿足利普希茨利普希茨(Lipschitz)條件條件，L稱為稱為y的的利普希茨常數(shù)利普希茨常數(shù)(簡(jiǎn)稱簡(jiǎn)稱Lips.常數(shù)常數(shù)).如果存在實(shí)數(shù)如果存在實(shí)數(shù)L0，使得，使得)3 . 1(.,),(),(212121RyyyyLyxfyxf 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 定理定理1 設(shè)設(shè)f在區(qū)域在區(qū)域D=(x, y)|a x b,y R上連續(xù)上連續(xù), 關(guān)關(guān)于于y滿足利普希茨條件，則對(duì)

4、任意滿足利普希茨條件，則對(duì)任意x0 a, b, y0 R，常，常微分方程初值問題微分方程初值問題(1.1)式和式和(1.2)式當(dāng)式當(dāng)x a, b時(shí)存在時(shí)存在唯一的連續(xù)可微解唯一的連續(xù)可微解y(x) .解的存在唯一性定理是常微分方程理論的基本內(nèi)解的存在唯一性定理是常微分方程理論的基本內(nèi)容，也是數(shù)值方法的出發(fā)點(diǎn)，此外還要考慮方程的解容，也是數(shù)值方法的出發(fā)點(diǎn)，此外還要考慮方程的解對(duì)擾動(dòng)的敏感性，它有以下結(jié)論對(duì)擾動(dòng)的敏感性，它有以下結(jié)論. 定理定理2 設(shè)設(shè)f在區(qū)域在區(qū)域D (如定理如定理1所定義所定義) 上連續(xù)上連續(xù), 且且關(guān)于關(guān)于y滿足利普希茨條件，設(shè)初值問題滿足利普希茨條件，設(shè)初值問題.)(),(

5、)(0sxyyxfxy 的解為的解為y(x, s)，則，則.),(),(21210ssesxysxyxxL 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)這個(gè)定理表明解對(duì)初值依賴的敏感性，它與右端這個(gè)定理表明解對(duì)初值依賴的敏感性，它與右端函數(shù)函數(shù)f有關(guān)，當(dāng)有關(guān)，當(dāng)f的的Lips.常數(shù)常數(shù)L比較小時(shí)，解對(duì)初值和比較小時(shí)，解對(duì)初值和右端函數(shù)相對(duì)不敏感，可視為好條件右端函數(shù)相對(duì)不敏感，可視為好條件. 若若L較大則可較大則可認(rèn)為壞條件，即病態(tài)問題認(rèn)為壞條件，即病態(tài)問題.如果右端函數(shù)可導(dǎo)，由中值定理有如果右端函數(shù)可導(dǎo)，由中值定理有.,),(),(),(212121之間之間在在yyyyyxfyxfyxf

6、 若假定若假定 在域在域D內(nèi)有界內(nèi)有界, 設(shè)設(shè) , 則則yyxf ),(Lyyxf ),(.),(),(2121yyLyxfyxf 它表明它表明f滿足利普希茨條件，且滿足利普希茨條件，且L的大小反映了右端函的大小反映了右端函數(shù)數(shù)f關(guān)于關(guān)于y變化的快慢，刻畫了初值問變化的快慢，刻畫了初值問(1.1)式和式和(1.2)式式是否為好條件是否為好條件. 這在數(shù)值求解中也是很重要的這在數(shù)值求解中也是很重要的.上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 雖然求解常微分方程有各種各樣的解析方法，但雖然求解常微分方程有各種各樣的解析方法，但解析方法只能用來(lái)求解一些特殊類型的方程，實(shí)際問解析方法只能用來(lái)求

7、解一些特殊類型的方程，實(shí)際問題中歸結(jié)出來(lái)的微分方程主要靠數(shù)值解法題中歸結(jié)出來(lái)的微分方程主要靠數(shù)值解法. 所謂所謂數(shù)值解法數(shù)值解法, 就是尋求解就是尋求解y(x)在一系列離散節(jié)點(diǎn)在一系列離散節(jié)點(diǎn) 121nnxxxx上的近似值上的近似值 y1,y2,yn,yn+1,. 相鄰兩個(gè)節(jié)點(diǎn)的間距相鄰兩個(gè)節(jié)點(diǎn)的間距hn=xn+1- -xn稱為稱為步長(zhǎng)步長(zhǎng). 今后如不特別說(shuō)明，總是假定今后如不特別說(shuō)明，總是假定 hi=h(i=1,2,)為為常數(shù)常數(shù), 這時(shí)節(jié)點(diǎn)為這時(shí)節(jié)點(diǎn)為xn=x0+nh(i=0,1,2,) (等距節(jié)點(diǎn)等距節(jié)點(diǎn)).上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 初值問題的初值問題的數(shù)值解法數(shù)

8、值解法有個(gè)有個(gè)基本特點(diǎn)基本特點(diǎn)，他們都采，他們都采取取“步進(jìn)式步進(jìn)式”，即求解過程順著節(jié)點(diǎn)排列的次序一步，即求解過程順著節(jié)點(diǎn)排列的次序一步一步地向前推進(jìn)一步地向前推進(jìn). 描述這類算法，只要給出用已知信描述這類算法，只要給出用已知信息息yn,yn- -1,yn- -2,計(jì)算計(jì)算yn+1的遞推公式的遞推公式. 本章首先要對(duì)常微分方程本章首先要對(duì)常微分方程(1.1)離散化，建立求解離散化，建立求解數(shù)值解的遞推公式數(shù)值解的遞推公式. 一類是計(jì)算一類是計(jì)算yn+1時(shí)只用到前一點(diǎn)的時(shí)只用到前一點(diǎn)的值值yn，稱為，稱為單步法單步法. 另一類是用到另一類是用到y(tǒng)n+1前面前面 k 點(diǎn)的值點(diǎn)的值yn,yn-1,

9、 yn-k+1，稱為，稱為k步法步法. 其次，要研究公式的其次，要研究公式的局局部截?cái)嗾`差部截?cái)嗾`差和和階階，數(shù)值解，數(shù)值解yn與精確解與精確解y(xn)的的誤差估計(jì)誤差估計(jì)及及收斂性收斂性，還有遞推公式的，還有遞推公式的計(jì)算穩(wěn)定性計(jì)算穩(wěn)定性等問題等問題.上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)9.2 簡(jiǎn)單的數(shù)值方法簡(jiǎn)單的數(shù)值方法9.2.1 歐拉法與后退歐拉法歐拉法與后退歐拉法 我們知道，在我們知道，在xy平面上，微分方程平面上，微分方程(1.1)式式的解的解y=f(x)稱作它的稱作它的積分曲線積分曲線，積分曲線積分曲線上一點(diǎn)上一點(diǎn)(x, y)的切的切線斜率等于函數(shù)線斜率等于函數(shù)f(x

10、, y)的值的值. 如果按如果按f(x, y)在在xy平面上平面上建立一個(gè)方向場(chǎng)，那么，建立一個(gè)方向場(chǎng)，那么，積分曲線積分曲線上每一點(diǎn)的切線上每一點(diǎn)的切線方向均與方向場(chǎng)在該點(diǎn)的方向相一致方向均與方向場(chǎng)在該點(diǎn)的方向相一致.基于上述幾何解釋，我們從初始點(diǎn)基于上述幾何解釋，我們從初始點(diǎn)P0(x0, y0)出發(fā)出發(fā),先依方向場(chǎng)在該點(diǎn)的方向推進(jìn)到先依方向場(chǎng)在該點(diǎn)的方向推進(jìn)到x=x1上一點(diǎn)上一點(diǎn)P1，然后，然后再?gòu)脑購(gòu)腜1點(diǎn)依方向場(chǎng)在該點(diǎn)的方向推進(jìn)到點(diǎn)依方向場(chǎng)在該點(diǎn)的方向推進(jìn)到 x=x2 上一點(diǎn)上一點(diǎn)P2 , 循環(huán)前進(jìn)做出一條循環(huán)前進(jìn)做出一條折線折線P0 P1 P2.上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下

11、頁(yè)下頁(yè)下頁(yè) 一般地，設(shè)已做出該折線的頂點(diǎn)一般地，設(shè)已做出該折線的頂點(diǎn)Pn, ,過過Pn(xn, yn)依方向場(chǎng)的方向再推進(jìn)到依方向場(chǎng)的方向再推進(jìn)到Pn+1(xn+1, yn+1)，顯然兩個(gè)，顯然兩個(gè)頂點(diǎn)頂點(diǎn)Pn, ,Pn+1的坐標(biāo)有關(guān)系的坐標(biāo)有關(guān)系),(,()(),(111nnnnnnnnnnnxyxfxyyxfhyyxxyy ) 1 .2 (),(1nnnnyxhfyy 這就是著名的這就是著名的( (顯式顯式) )歐拉歐拉( (Euler) )公式公式. . 若初值若初值y0已已知，則依公式知，則依公式(2.1)可逐次逐步算出各點(diǎn)數(shù)值解可逐次逐步算出各點(diǎn)數(shù)值解. .即即),(0001yxhf

12、yy ),(1112yxhfyy 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 例例1 用歐拉公式求解初值問題用歐拉公式求解初值問題)2 . 2(. 1)0(),10(2 yxyxyy 解解 取步長(zhǎng)取步長(zhǎng)h=0.1，歐拉公式的具體形式為，歐拉公式的具體形式為)2(1nnnnnyxyhyy 其中其中xn=nh=0.1n (n=0,1,10), 已知已知y0 =1, 由此式可得由此式可得191818. 1)1 . 12 . 01 . 1 ( 1 . 01 . 1)2(1 . 11 . 01)2(1111200001 yxyhyyyxyhyy上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)依次計(jì)

13、算下去，依次計(jì)算下去，部分計(jì)算結(jié)果部分計(jì)算結(jié)果見下表見下表. xy21 與準(zhǔn)確解與準(zhǔn)確解 相比，可看出歐拉公式的計(jì)算結(jié)相比，可看出歐拉公式的計(jì)算結(jié)果精度很差果精度很差. xn 歐拉公式數(shù)值解歐拉公式數(shù)值解yn準(zhǔn)確解準(zhǔn)確解y(xn) 誤差誤差 0.2 0.4 0.6 0.8 1.0 1.191818 1.358213 1.508966 1.649783 1.784770 1.183216 1.341641 1.483240 1.612452 1.732051 0.008602 0.016572 0.025726 0.037331 0.052719上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)

14、歐拉公式具有明顯的幾何意義歐拉公式具有明顯的幾何意義, , 就是就是用折線近似用折線近似代替方程的解曲線代替方程的解曲線，因而常稱公式，因而常稱公式(2.1)為為歐拉折線法歐拉折線法. .( )yy xxynx1nxnp1np1np x 還可以通過幾何直觀來(lái)考察歐拉方法的精度還可以通過幾何直觀來(lái)考察歐拉方法的精度. .假假設(shè)設(shè)yn=y(xn), ,即頂點(diǎn)即頂點(diǎn)Pn落在積分曲線落在積分曲線y=y(x)上，那么，上，那么，按歐拉方法做出的折線按歐拉方法做出的折線PnPn+1便是便是y=y(x)過點(diǎn)過點(diǎn)Pn的切線的切線. .從圖形上看從圖形上看, ,這這樣定出的頂點(diǎn)樣定出的頂點(diǎn)Pn+1顯著顯著地偏離

15、了原來(lái)的積分曲地偏離了原來(lái)的積分曲線，可見歐拉方法是線，可見歐拉方法是相相當(dāng)粗糙當(dāng)粗糙的的. .上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 為了分析計(jì)算公式的精度，通?？捎锰├照归_為了分析計(jì)算公式的精度，通?？捎锰├照归_將將y(xn+1)在在xn處展開，則有處展開，則有).,()(2),()()(2)()()()(1221 nnnnnnnnnnnnxxyhyxhfxyyhxyhxyhxyxy 在在yn=y(xn)的前提下，的前提下，f(xn,yn )=f(xn,y(xn)=y ( (xn n) ). .于是于是可得歐拉法可得歐拉法(2.1)的的公式誤差公式誤差為為) 3 . 2(),(2

16、)(2)(2211nnnnxyhyhyxy 稱為此方法的稱為此方法的局部截?cái)嗾`差局部截?cái)嗾`差. .上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 如果對(duì)方程如果對(duì)方程(1.1)從從xn到到xn+1積分，得積分，得) 4 . 2(.)(,()()(11 nnxxnndttytfxyxy右端積分用右端積分用左矩形公式左矩形公式hf(xn,y(xn)近似，再以近似，再以yn代替代替y(xn)，yn+1代替代替y(xn+1)也得到歐拉公式也得到歐拉公式(2.1)，局部截，局部截?cái)嗾`差也是斷誤差也是(2.3).稱為稱為( (隱式隱式) )后退的歐拉公式后退的歐拉公式. .它也可以通過利用均差它也可以通

17、過利用均差近似導(dǎo)數(shù)近似導(dǎo)數(shù)y (xn+1)，即，即 如果右端積分用如果右端積分用右矩形公式右矩形公式hf(xn+1,y(xn+1)近似，近似，則得到另一個(gè)公式則得到另一個(gè)公式) 5 . 2 (),(111 nnnnyxhfyy直接得到直接得到. .),(,()()()(1111 nnnnnnnxyxfxyxxxyxy上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 后退的歐拉公式與歐拉公式有著本質(zhì)的區(qū)別后退的歐拉公式與歐拉公式有著本質(zhì)的區(qū)別, 后后者是關(guān)于者是關(guān)于yn+1的一個(gè)直接計(jì)算公式，這類公式稱作是的一個(gè)直接計(jì)算公式，這類公式稱作是顯式的顯式的；前者公式的；前者公式的右端含有未知的右端含

18、有未知的yn+1，它實(shí)際上，它實(shí)際上是關(guān)于是關(guān)于yn+1的一個(gè)函數(shù)方程的一個(gè)函數(shù)方程, ,這類方程稱作是這類方程稱作是隱式的隱式的. . 顯式顯式與與隱式隱式兩類方法各有特點(diǎn)，考慮到數(shù)值穩(wěn)兩類方法各有特點(diǎn)，考慮到數(shù)值穩(wěn)定性等其他因素，人們有時(shí)需要選用定性等其他因素，人們有時(shí)需要選用隱式隱式方法，但方法，但使用使用顯式顯式算法遠(yuǎn)比算法遠(yuǎn)比隱式隱式方便方便. . 隱式方程隱式方程(2.5)通常用迭代法求解，而迭代過程通常用迭代法求解，而迭代過程的實(shí)質(zhì)是的實(shí)質(zhì)是逐步逐步顯式化顯式化. .上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 設(shè)用歐拉公式設(shè)用歐拉公式),() 0(1nnnnyxhfyy 給

19、出迭代初值給出迭代初值 ，用它代入，用它代入(2.5)式的式的右端，使之轉(zhuǎn)右端，使之轉(zhuǎn)化為顯式，直接計(jì)算得化為顯式，直接計(jì)算得) 0(1 ny),() 0(11) 1 (1 nnnnyxhfyy然后再用然后再用 代入代入(2.5)式，又有式，又有) 1 (1 ny).,() 1 (11) 2(1 nnnnyxhfyy如此反復(fù)進(jìn)行，得如此反復(fù)進(jìn)行，得) 6 . 2 ()., 1 , 0(),()(11) 1(1 kyxhfyyknnnkn上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)由于由于f(x, y)對(duì)對(duì)y滿足滿足Lipschitz條件條件(1.3). 由由(2.6)減減(2.5)得得),

20、(),(11)(111) 1(1 nnknnnknyxfyxfhyy.1)(1 nknyyhL由此可知，只要由此可知，只要hL1，迭代法，迭代法(2.6)就收斂到解就收斂到解yn+1. .關(guān)于后退歐拉方法的公式誤差，從積分公式看到它關(guān)于后退歐拉方法的公式誤差，從積分公式看到它與歐拉法是相似的與歐拉法是相似的. .上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)9.2.2 梯形方法梯形方法 為得到比為得到比歐拉法精度高的計(jì)算公式，在等歐拉法精度高的計(jì)算公式，在等式式(2.4) 右端積分用梯形求積公式近似右端積分用梯形求積公式近似, 并用并用yn代替代替y(xn), yn+1代替代替y(xn+1)

21、，則得，則得 ) 7 . 2(,),(),(2111 nnnnnnyxfyxfhyy稱為稱為矩形方法矩形方法. 矩形方法矩形方法是是隱式單步法隱式單步法，用迭代法求解，同后，用迭代法求解，同后退的歐拉方法一樣，仍用歐拉法提供迭代初值，則退的歐拉方法一樣，仍用歐拉法提供迭代初值，則矩形迭代公式矩形迭代公式為為上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) ) 8 . 2()., 1 , 0(),(),(2);,()(11) 1()0(11 kyxfyxfhyyyxhfyyknnnnnknnnnn為了分析迭代過程的收斂性為了分析迭代過程的收斂性, 將將(2.7)與與(2.8)相減相減, 得得),

22、(),(2)(1111) 1(11knnnnknnyxfyxfhyy 于是有于是有,2)(11) 1(11knnknnyyhLyy 使得使得, 12 hL則當(dāng)則當(dāng)k時(shí)有時(shí)有 , 這說(shuō)明迭代過程這說(shuō)明迭代過程(2.8)是收斂的是收斂的.1)(1 nknyy上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)9.2.3 改進(jìn)的歐拉公式改進(jìn)的歐拉公式我們看到，梯形方法雖然提高了精度，但其算法我們看到，梯形方法雖然提高了精度，但其算法復(fù)雜，在應(yīng)用迭代公式復(fù)雜，在應(yīng)用迭代公式(2.8)進(jìn)行實(shí)際計(jì)算時(shí)，每迭代進(jìn)行實(shí)際計(jì)算時(shí)，每迭代一次，都要重新計(jì)算函數(shù)一次，都要重新計(jì)算函數(shù)f(x, y )的值，而迭代又要反的

23、值，而迭代又要反復(fù)進(jìn)行若干次，計(jì)算量很大，而且往往難以預(yù)測(cè)復(fù)進(jìn)行若干次，計(jì)算量很大，而且往往難以預(yù)測(cè). 為為了控制計(jì)算量，通常只迭代一兩次就轉(zhuǎn)入下一步的計(jì)了控制計(jì)算量，通常只迭代一兩次就轉(zhuǎn)入下一步的計(jì)算，這算，這就簡(jiǎn)化了算法就簡(jiǎn)化了算法. .具體地說(shuō)，我們先用歐拉公式求得一個(gè)初步的近具體地說(shuō)，我們先用歐拉公式求得一個(gè)初步的近似值似值 ，稱之為，稱之為預(yù)測(cè)值預(yù)測(cè)值，此預(yù)測(cè)值，此預(yù)測(cè)值 的精度可能很的精度可能很差，再用梯形公式差，再用梯形公式(2.7)將它校正一次，即按將它校正一次，即按(2.8)式迭式迭代一次得代一次得yn+1，這個(gè)結(jié)果，這個(gè)結(jié)果稱之為稱之為校正值校正值. .1 ny1 ny上頁(yè)

24、上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)這樣建立的這樣建立的預(yù)測(cè)預(yù)測(cè)- -校正系統(tǒng)校正系統(tǒng)通常稱為通常稱為改進(jìn)的歐拉公式改進(jìn)的歐拉公式:或表示為下列或表示為下列平均化形式平均化形式 ).(21),(),(11cpnpnncnnnpyyyyxhfyyyxhfyy( (2.9) ), ),(1nnnnyxhfyy 預(yù)測(cè)預(yù)測(cè).),(),(2111 nnnnnnyxfyxfhyy校正校正上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 例例2 用改進(jìn)的歐拉法用改進(jìn)的歐拉法解解例例1中的初值問題中的初值問題(2.2). . 解解 仍取步長(zhǎng)仍取步長(zhǎng)h=0.1, ,改進(jìn)的歐拉公式為改進(jìn)的歐拉公式為 ).

25、(21),2(),2(11cpnpnpncnnnnpyyyyxyhyyyxyhyy部分計(jì)算結(jié)果部分計(jì)算結(jié)果見下表見下表 xnyn 誤差誤差 xnyn 誤差誤差00.20.4 1. 1.184096 1.34336000.0000880.0017190.60.81.01.4859561.6164761.7378690.0027160.0040240.05818同例同例1中的歐拉法的計(jì)算結(jié)果比較，明顯改善了精度中的歐拉法的計(jì)算結(jié)果比較，明顯改善了精度.上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)9.2.4 單步法的局部截?cái)嗾`差與階單步法的局部截?cái)嗾`差與階 初值問題初值問題(1.1),(1.2)

26、的的單步法單步法可用可用一般形式一般形式表示為表示為)10. 2(, ),(11hyyxhyynnnnn 其中多元函數(shù)其中多元函數(shù) 與與f(x, y )有關(guān)，當(dāng)有關(guān)，當(dāng) 含有含有yn+1時(shí)，方法時(shí)，方法是是隱式隱式的，若不含的，若不含yn+1則為則為顯式方法顯式方法，所以顯式單步，所以顯式單步法可表示為法可表示為 (x, y, h)稱為稱為增量函數(shù)增量函數(shù)，例如對(duì)歐拉法，例如對(duì)歐拉法(2.1)有有)11. 2(, ),(1hyxhyynnnn . ),(),(yxfhyxnn 它的局部截?cái)嗾`差已由它的局部截?cái)嗾`差已由(2.3)給出給出, 對(duì)一般顯式單步法對(duì)一般顯式單步法則可如下定義則可如下定義

27、. .上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 定義定義1 設(shè)設(shè)y(x)是初值問題是初值問題(1.1),(1.2)的的準(zhǔn)確解準(zhǔn)確解, 稱稱)12. 2().),(,(11hxyxhyyTnnnnn 為顯式單步法為顯式單步法( (2.11) )的的局部截?cái)嗾`差局部截?cái)嗾`差. . Tn+1之所以稱為之所以稱為局部的局部的，是假設(shè)在，是假設(shè)在xn前各步?jīng)]有前各步?jīng)]有誤差誤差. .當(dāng)當(dāng)yn= =y(xn)時(shí)，計(jì)算一步，則有時(shí)，計(jì)算一步，則有.),(,()()(),()()(11111 nnnnnnnnnnnThxyxhxyxyhyxhyxyyxy 所以，局部截?cái)嗾`差可理解為用方法所以，局部截?cái)?/p>

28、誤差可理解為用方法(2.11)計(jì)算計(jì)算一步的誤差，也即公式一步的誤差，也即公式(2.11)中用準(zhǔn)確解中用準(zhǔn)確解y(x)代替數(shù)代替數(shù)值解產(chǎn)生的公式誤差值解產(chǎn)生的公式誤差. 根據(jù)定義根據(jù)定義, 顯然顯然歐拉法的局部歐拉法的局部截?cái)嗾`差截?cái)嗾`差為為上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)).()(2)()()()(,()()(3211hOxyhxyhxyhxyxyxhfxyxyTnnnnnnnnn 即為即為(2.3)的結(jié)果的結(jié)果. 這里這里 稱為稱為局部截?cái)嗾`差主局部截?cái)嗾`差主項(xiàng)項(xiàng). 顯然顯然Tn+1= =O(h2). 一般情形的定義如下一般情形的定義如下)(22nxyh 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上

29、頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 定義定義2 設(shè)設(shè)y (x)是初值問題的準(zhǔn)確解，若存在最是初值問題的準(zhǔn)確解，若存在最大整數(shù)大整數(shù)p使顯式單步法使顯式單步法( (2.11) )的的局部截?cái)嗾`差局部截?cái)嗾`差滿足滿足)13. 2().(),()()(11 pnhOhyxhxyhxyT 則稱方法則稱方法( (2.11) )具有具有p階精度階精度. .若將若將( (2.11) )展開式寫成展開式寫成).()(,(211 ppnnnhOhxyxT 則則 稱為稱為局部截?cái)嗾`差主項(xiàng)局部截?cái)嗾`差主項(xiàng). 1)(,( pnnhxyx 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)).()(2)()()()()(2)

30、()(,()()(322321111hOxyhhOxyhxyhhOxyhxyhxyxhfxyxyTnnnnnnnnnn 以上定義對(duì)隱式單步法以上定義對(duì)隱式單步法(2.10)也是適用的也是適用的. .例如，例如，對(duì)后退歐拉法對(duì)后退歐拉法(2.5)其局部截?cái)嗾`差為其局部截?cái)嗾`差為這里這里p=1是是1階方法階方法，局部截?cái)嗾`差主項(xiàng)局部截?cái)嗾`差主項(xiàng)為為).(22nxyh 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)同樣對(duì)梯形法同樣對(duì)梯形法( (2.7) )有有).()(12)()(2)()()(2)(! 3)(2)()()(2)()(434232111hOxyhhOxyhxyhxyxyhxyhxy

31、hxyhxyxyhxyxyTnnnnnnnnnnnnn 所以梯形方法所以梯形方法(2.7)是是二階的二階的. 其其局部截?cái)嗾`差主項(xiàng)局部截?cái)嗾`差主項(xiàng)為為).(124nxyh 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)9.3 龍格龍格- -庫(kù)塔方法庫(kù)塔方法 對(duì)許多實(shí)際問題來(lái)說(shuō)，歐拉公式與改進(jìn)歐拉對(duì)許多實(shí)際問題來(lái)說(shuō)，歐拉公式與改進(jìn)歐拉公式精度還不能滿足要求，為此從另一個(gè)角度來(lái)分公式精度還不能滿足要求，為此從另一個(gè)角度來(lái)分析這兩個(gè)公式的特點(diǎn)，從而探索一條構(gòu)造高精度方析這兩個(gè)公式的特點(diǎn)，從而探索一條構(gòu)造高精度方法的途徑法的途徑. 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)9.3.1 顯式龍格

32、顯式龍格- -庫(kù)塔法的一般形式庫(kù)塔法的一般形式 上節(jié)給出了顯式單步法的表達(dá)式上節(jié)給出了顯式單步法的表達(dá)式( (2.11) ), 其局部其局部截?cái)嗾`差為截?cái)嗾`差為O(hp+1)，對(duì)歐拉法，對(duì)歐拉法Tn+1=O(h2)，即方法為，即方法為p=1階階，若用改進(jìn)的歐拉法，若用改進(jìn)的歐拉法( (2.9) )式，它可表示為式，它可表示為) 1 . 3().,(,(),(21nnnnnnnnyxhfyhxfyxfhyy 此時(shí)增量函數(shù)為此時(shí)增量函數(shù)為) 2 . 3().,(,(),(21),(nnnnnnnnyxhfyhxfyxfhyx 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)它比歐拉法的它比歐拉法的

33、 (xn, yn, h)= =f(xn, yn), 增加了計(jì)算一個(gè)右增加了計(jì)算一個(gè)右函數(shù)函數(shù)f 的值，可望的值，可望 p=2. .若要使得到的公式階數(shù)若要使得到的公式階數(shù)p更大更大, 就必須包含更多的就必須包含更多的f 值值. . 實(shí)際上從方程實(shí)際上從方程(1.1)等價(jià)的等價(jià)的積分形式積分形式(2.4) ，即，即)3 . 3(,)(,()()(11 nnxxnndxxyxfxyxy若要使公式階數(shù)提高，就必須使右端積分的數(shù)值求若要使公式階數(shù)提高，就必須使右端積分的數(shù)值求積公式精度提高，它必然要增加求積節(jié)點(diǎn)，為此可積公式精度提高，它必然要增加求積節(jié)點(diǎn)，為此可將將( (3.3) )的右端用求積公式表

34、示為的右端用求積公式表示為.)(,()(,(11 riininixxhxyhxfchdxxyxfnn 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)一般說(shuō)來(lái)，點(diǎn)數(shù)一般說(shuō)來(lái)，點(diǎn)數(shù)r 越多，精度越高，上式右端相當(dāng)于越多，精度越高，上式右端相當(dāng)于增量函數(shù)增量函數(shù) (x, y, h)，為得到便于計(jì)算的顯式方法，為得到便于計(jì)算的顯式方法, 可可類似于改進(jìn)歐拉法類似于改進(jìn)歐拉法(3.1),(3.2)，將公式表示為將公式表示為) 4 . 3(, ),(1hyxhyynnnn 其中其中) 5 . 3(,),(1 riiinnKchyx , ),(1nnyxfK , 2).,(11riKhyhxfKijjij

35、nini 這里這里ci, i, ij均為常數(shù)均為常數(shù). (3.4)和和(3.5)稱為稱為r級(jí)顯式級(jí)顯式龍格龍格- -庫(kù)塔庫(kù)塔(Runge- -Kutta)法法, 簡(jiǎn)稱簡(jiǎn)稱R- -K方法方法.上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 當(dāng)當(dāng)r=1, (xn, yn, h)= =f(xn, yn)時(shí)，就是歐拉法，此時(shí)，就是歐拉法，此時(shí)方法的階為時(shí)方法的階為p=1. 當(dāng)當(dāng)r=2時(shí)，改進(jìn)歐拉法時(shí)，改進(jìn)歐拉法(3.1)式就是式就是其中的一種其中的一種, 下面將證明其階下面將證明其階p=2. 要使公式要使公式(3.4),(3.5)具有更高的階具有更高的階p，就要增加點(diǎn)數(shù)，就要增加點(diǎn)數(shù)r. 下面我們只

36、就下面我們只就r=2推導(dǎo)推導(dǎo)R- -K方法方法. 并給出并給出 r=3,4 時(shí)的常用公式，其推導(dǎo)時(shí)的常用公式，其推導(dǎo)方法與方法與r=2時(shí)類似，只是推導(dǎo)計(jì)算較復(fù)雜時(shí)類似，只是推導(dǎo)計(jì)算較復(fù)雜.上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)9.3.2 二階顯式二階顯式R- -K方法方法 對(duì)對(duì)r=2的的R- -K方法方法，由，由(3.4),(3.5)式可得式可得如下如下計(jì)算計(jì)算公式公式)6 . 3().,(),(),(12122122111 hKyhxfKyxfKKcKchyynnnnnn 這里這里 c1, c2, 2, 21 均為待定常數(shù)，我們希望適當(dāng)選取均為待定常數(shù)，我們希望適當(dāng)選取這些系數(shù)，使

37、公式階數(shù)這些系數(shù)，使公式階數(shù) p 盡量高盡量高. 根據(jù)局部截?cái)嗾`差根據(jù)局部截?cái)嗾`差定義，推導(dǎo)出定義，推導(dǎo)出(3.6)的局部截?cái)嗾`差為的局部截?cái)嗾`差為上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)) 7 . 3(),(),()()(2122111nnnnnnnnhfyhxfcyxfchxyxyT ),(! 32)(4321hOyhyhyhyxynnnnn 其中其中這里這里yn= =y(xn), yn+1= =y(xn+1). 為得到為得到Tn+1的階的階p，要將上，要將上式各項(xiàng)在式各項(xiàng)在(xn, yn)處做泰勒展開，由于處做泰勒展開，由于f(x, y )是二元函是二元函數(shù)，故要用二元泰勒展開，

38、各項(xiàng)展開式為數(shù)，故要用二元泰勒展開，各項(xiàng)展開式為)8 . 3(.2),(),()(,(,),(2 ynxyyynxynxxnnnynnnxnnnnnnnfffffffffyyxffyxfxyxfdxdyfyxfy上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)).(),(),(),(2212212hOhfyxfhyxffhfyhxfnnnynnxnnnn 將以上結(jié)果代入將以上結(jié)果代入(3.7)，則有，則有)()()(221111KcKchxyxyTnnn 22221),()21()1 (hyxfchfccnnxn ).(),()21(32212hOhfyxfcnnny ),(),()()(21

39、22111nnnnnnnnhfyhxfcyxfchxyxyT )()()(232122132hOhffhffcfchhOfffhhfnyxnnnyxn 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)要使公式要使公式(3.6)具有具有p=2階，必須使階，必須使12222211110,0,0.(3.9)22cccc 222211211,1.22cccc 即即(3.9)的的解是不唯一的解是不唯一的. . 可令可令c2=a0，則得，則得122111,.2caa 這樣得到的公式這樣得到的公式稱為稱為二階二階R- -K方法方法.上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 則由此可以看出在改進(jìn)的歐拉

40、公式中相當(dāng)于取則由此可以看出在改進(jìn)的歐拉公式中相當(dāng)于取(xn,yn), (xn+1,yn+1)兩兩點(diǎn)處點(diǎn)處斜率的平均值斜率的平均值，近似代替平，近似代替平均斜率，其精度比歐拉公式提高了均斜率，其精度比歐拉公式提高了. . ).,(),(),(21121211hKyxfKyxfKKKhyynnnnnn 如取如取a=1/2，則，則c1= c2=1/2, 2= 21=1. . 這就是這就是改進(jìn)改進(jìn)的歐拉公式的歐拉公式(3.1)式式.上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)稱為稱為中點(diǎn)公式中點(diǎn)公式( (變形的歐拉公式變形的歐拉公式) )，相當(dāng)于數(shù)值積分的，相當(dāng)于數(shù)值積分的中矩形公式中矩形公式.

41、 .也可以表示為也可以表示為)10. 3().2,2(),(,12121 KhyhxfKyxfKhKyynnnnnn 如取如取a=1，則，則c1=0, c2=1, 2= 21=1/2. . 得計(jì)算公式得計(jì)算公式.),(2,2(1nnnnnnyxfhyhxhfyy 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 對(duì)對(duì)r=2的的R- -K公式公式(3.6)能否使局部誤差提高到能否使局部誤差提高到O(h4)? 為此為此 需把需把K2多展開一項(xiàng)，從多展開一項(xiàng)，從(3.8)的的 看到展看到展開式中的項(xiàng)開式中的項(xiàng) 是不能通過選擇參數(shù)消掉是不能通過選擇參數(shù)消掉的，實(shí)際上要使的，實(shí)際上要使 h3 的項(xiàng)為零，

42、需增加的項(xiàng)為零，需增加3個(gè)方程，要確個(gè)方程，要確定定4個(gè)參數(shù)個(gè)參數(shù)c1, c2, 2及及 21，這是不可能的，這是不可能的. 故故r=2的顯的顯式式R- -K方法的階只能是方法的階只能是p=2，而不能得到三階公式，而不能得到三階公式.ffffyxy2)( ny 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)9.3.3 三階與四階顯式三階與四階顯式R- -K方法方法 要得到三階顯式要得到三階顯式R- -K方法方法，必須，必須r=3. 此時(shí)計(jì)算此時(shí)計(jì)算(3.4), (3.5)的公式的公式表示為表示為)11. 3().,(),(),(),(232131331212213322111 hKhKyhx

43、fKhKyhxfKyxfKKcKcKchyynnnnnnnn 其中其中c1, c2, c3及及 2, 21, 3, 31, 32均為待定常數(shù)，公式均為待定常數(shù)，公式(3.11)的局部截?cái)嗾`差為的局部截?cái)嗾`差為)()()(33221111KcKcKchxyxyTnnn 上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)只要只要K1, K2將按二元泰勒展開，使將按二元泰勒展開，使Tn+1O(h4)，可得，可得待定參數(shù)滿足方程待定參數(shù)滿足方程)12. 3(.61,31,21, 13223233222332232313212321 cccccccc上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)這是這是

44、8個(gè)未知數(shù)個(gè)未知數(shù)6個(gè)方程的方程組，解不是唯一的個(gè)方程的方程組，解不是唯一的. 可可以得到很多公式以得到很多公式. 滿足條件滿足條件( (3.12) )的公式的公式( (3.11) )統(tǒng)稱統(tǒng)稱為為三階三階R- -K公式公式. 下面只給出下面只給出其中一個(gè)常見的公式其中一個(gè)常見的公式. ).2,(),2,2(),(),4(62131213211hKhKyhxfKKhyhxfKyxfKKKKhyynnnnnnnn此公式稱為此公式稱為庫(kù)塔三階方法庫(kù)塔三階方法. .上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 繼續(xù)上述過程繼續(xù)上述過程, 經(jīng)過較復(fù)雜的數(shù)學(xué)演算經(jīng)過較復(fù)雜的數(shù)學(xué)演算, 可以導(dǎo)可以導(dǎo)出出

45、各種四階各種四階R- -K公式公式,下列下列經(jīng)典公式經(jīng)典公式是其中常用的一是其中常用的一個(gè)：個(gè)： ).,(),2,2(),2,2(),(),22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn 四階四階R- -K方法方法的每一步需要計(jì)算四次函數(shù)值的每一步需要計(jì)算四次函數(shù)值f，可以證明其局部截?cái)嗾`差為可以證明其局部截?cái)嗾`差為O(h5). . 證明從略證明從略. .上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 例例3 設(shè)取步長(zhǎng)設(shè)取步長(zhǎng)h=0.2, ,從從x=0直到直到x=1用四階龍格用四階龍格- -庫(kù)塔方法求解庫(kù)塔方法求解例例1中

46、的初值問題中的初值問題(2.2)式式. . 解解 這里這里, ,經(jīng)典的四階龍格經(jīng)典的四階龍格- -庫(kù)塔公式庫(kù)塔公式(3.13)具有具有 .)(2,22,22,2),22(633422231212143211hKyhxhKyKKyhxKhyKKyhxKhyKyxyKKKKKhyynnnhnnnhnnnnnnnn上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)計(jì)算結(jié)果計(jì)算結(jié)果見下表見下表 xnyny(xn) xnyny(xn) 00.20.4 1. 1.1832 1.3417 0 1.1832 3.24160.60.81.0 1.4833 1.6125 1.7321 1.4832 1.6125 1

47、.7321比較例比較例3和例和例2的計(jì)算結(jié)果，顯然以龍格的計(jì)算結(jié)果，顯然以龍格- -庫(kù)塔方庫(kù)塔方法的精度為高法的精度為高. 要注意，雖然四階龍格要注意，雖然四階龍格- -庫(kù)塔方法的庫(kù)塔方法的計(jì)算量計(jì)算量(每一步要每一步要4次計(jì)算函數(shù)次計(jì)算函數(shù)f )比改進(jìn)的歐拉方法比改進(jìn)的歐拉方法(它是一種二階龍格它是一種二階龍格- -庫(kù)塔方法，每一步只要庫(kù)塔方法，每一步只要2次計(jì)算次計(jì)算函數(shù)函數(shù)f )大一倍，但由于這里放大了步長(zhǎng)大一倍，但由于這里放大了步長(zhǎng)(h=0.2)，兩個(gè)，兩個(gè)例子的算法所耗費(fèi)的計(jì)算量幾乎相同例子的算法所耗費(fèi)的計(jì)算量幾乎相同. 這個(gè)例子又一這個(gè)例子又一次顯示了選擇算法的重要意義次顯示了選擇

48、算法的重要意義.上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè) 然而值得指出的是，龍格然而值得指出的是，龍格- -庫(kù)塔方法的推導(dǎo)基于庫(kù)塔方法的推導(dǎo)基于泰勒展開泰勒展開方法，因而它要求所求的方法，因而它要求所求的解解具有較好的光具有較好的光滑性質(zhì)滑性質(zhì). 反之，如果反之，如果解解的光滑性差，那么，使用龍的光滑性差，那么，使用龍格格- -庫(kù)塔方法求得的數(shù)值解，其精度可能反而不如改庫(kù)塔方法求得的數(shù)值解，其精度可能反而不如改進(jìn)的歐拉方法進(jìn)的歐拉方法. 實(shí)際計(jì)算時(shí)，我們應(yīng)當(dāng)實(shí)際計(jì)算時(shí)，我們應(yīng)當(dāng)針對(duì)問題的具針對(duì)問題的具體特點(diǎn)選擇合適的算法體特點(diǎn)選擇合適的算法.上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)

49、下頁(yè)微分方程組和高階方程初值問題的數(shù)值解微分方程組和高階方程初值問題的數(shù)值解 歐拉方法和龍格歐拉方法和龍格-庫(kù)塔方法可直接推廣到微分方程組庫(kù)塔方法可直接推廣到微分方程組0000( , , )( , , )(), ()yf x y zzg x y zy xyz xz向前歐拉公式向前歐拉公式 ),(yyxfy 1yy),(21221yyxfyyy高階方程需要先降階為一階微分方程組高階方程需要先降階為一階微分方程組 , 2 , 1 , 0),(),(11nzyxhgzzzyxhfyynnnnnnnnnn上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)上頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)下頁(yè)龍格龍格庫(kù)塔方法的庫(kù)塔方法的 MATLAB 實(shí)現(xiàn)實(shí)現(xiàn) TnTnfffxxxxtxxtftx),(,),(,)(),()(1100t,x=ode23(f,ts,x0,opt) 3級(jí)級(jí)2階龍格階龍格-庫(kù)塔公式庫(kù)塔公式 t,x=ode45(f,ts,x0,opt) 5級(jí)級(jí)4階龍格階龍格-庫(kù)塔公式庫(kù)塔公式 f是待解方程寫成是待解方程寫成的函數(shù)的函數(shù)m文件：文件： function dx=f(t,x)dx=f1; f2; fn;ts