第六章常微分方程初值問題的數(shù)值解法n

1、常微分方程初值問題常微分方程初值問題的數(shù)值解法的數(shù)值解法第第6章章引言引言在實(shí)際問題中，常需要求解微分方程在實(shí)際問題中，常需要求解微分方程(如發(fā)電機(jī)運(yùn)動(dòng)方程如發(fā)電機(jī)運(yùn)動(dòng)方程)。只有簡(jiǎn)單的和典型的微分方程可以求出解析解，而在只有簡(jiǎn)單的和典型的微分方程可以求出解析解，而在實(shí)際問題中的微分方程往往無法求出解析解。實(shí)際問題中的微分方程往往無法求出解析解。在高等數(shù)學(xué)中我們見過以下常微分方程：在高等數(shù)學(xué)中我們見過以下常微分方程： 0)(),(yaybxayxfy-(1) )(,)(),(0ayyaybxayyxfy-(2)一階常微分方程一階常微分方程 nybyyaybxayyxfy)(,)(),(0-(3

2、)(1),(2)式稱為式稱為初值問題初值問題,(,(3)式稱為式稱為邊值問題邊值問題 2002212210012111)(),()(),(yxyyyxfyyxyyyxfy-(4)另外另外, ,在實(shí)際應(yīng)用中還經(jīng)常需要求解常微分方程組在實(shí)際應(yīng)用中還經(jīng)常需要求解常微分方程組: 本課程主要研究問題本課程主要研究問題一階常微分方程一階常微分方程(1)的數(shù)值解法的數(shù)值解法,我們首先介紹初值問題我們首先介紹初值問題(1)(1)的解存在的條件的解存在的條件定理定理 只要只要 f (x, y) 連續(xù)，且關(guān)于連續(xù)，且關(guān)于 y 滿足滿足 Lipschitz 條件條件，即存在與即存在與 x, y 無關(guān)的常數(shù)無關(guān)的常數(shù)

3、 L 使使對(duì)任意定義在對(duì)任意定義在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，則初值問都成立，則初值問題（題（1）存在唯一解存在唯一解。| ),(),(|2121yyLyxfyxf bxxxxan 210), 2 , 1()(nkyxykk 的近似值的近似值上函數(shù)值上函數(shù)值的數(shù)值解的數(shù)值解就是問題就是問題而而)1(), 2 , 1(nkyk 上的一系列離散點(diǎn)上的一系列離散點(diǎn)在區(qū)間在區(qū)間就是求未知函數(shù)就是求未知函數(shù),)(baxy（通常采用（通常采用等距節(jié)點(diǎn)）等距節(jié)點(diǎn)）對(duì)于問題對(duì)于問題(1)(1) 0)(),(yaybxayxfy要求它的要求它的數(shù)值解數(shù)值解常微分方程數(shù)值解公式的推

4、導(dǎo)常微分方程數(shù)值解公式的推導(dǎo) 求初值問題數(shù)值解的方法是求初值問題數(shù)值解的方法是步進(jìn)法步進(jìn)法，即從已知的初值，即從已知的初值y0出發(fā)，通過出發(fā)，通過一定的計(jì)算一定的計(jì)算求求y1 ，然后由，然后由y1或或y0和和y1求出求出y2 ，依次計(jì)算到依次計(jì)算到y(tǒng)n ，即在計(jì)算出，即在計(jì)算出yk后計(jì)算后計(jì)算yk+1 ，這時(shí)有，這時(shí)有單步法單步法：計(jì)算：計(jì)算yk+1時(shí)，只利用時(shí)，只利用yk多步法多步法：計(jì)算：計(jì)算yk+1時(shí)，用到時(shí)，用到y(tǒng)k, yk-1, yk-2,常微分方程數(shù)值解公式常微分方程數(shù)值解公式的主要推導(dǎo)方法的主要推導(dǎo)方法泰勒展開泰勒展開利用差商利用差商利用數(shù)值積分法利用數(shù)值積分法1、泰勒展開的求解

5、思路、泰勒展開的求解思路：將將 按泰勒級(jí)數(shù)展開按泰勒級(jí)數(shù)展開 hxyxykk 1 kkkkyhxhyxyxy 21! 2)( 211() , ()kkkkkkkhy xy xhy xy xy xhfxy x 略略去去 1, 1 , 0,10nkyxhfyyyaykkkk用用 的近似值的近似值 代入上式右端，記所得結(jié)果代入上式右端，記所得結(jié)果為為，則得到數(shù)值解序列的計(jì)算公式，則得到數(shù)值解序列的計(jì)算公式：()ky xky1ky 2、化導(dǎo)數(shù)為差商的求解方法思路：、化導(dǎo)數(shù)為差商的求解方法思路：若在點(diǎn)若在點(diǎn) 處的導(dǎo)數(shù)用差商來近似代替，如向前差商處的導(dǎo)數(shù)用差商來近似代替，如向前差商 hxyxyxykkk

6、1kx則微分方程初值問題化為則微分方程初值問題化為 1, 1 , 001nkyayxyhxyxykkk將近似號(hào)改為等號(hào)，精確解將近似號(hào)改為等號(hào)，精確解 改為近似解改為近似解 ，得，得ky kxy 1, 1 , 0,10nkyxhfyyyaykkkk3、數(shù)值積分的求解思路、數(shù)值積分的求解思路：如果將微分方程如果將微分方程 在各小區(qū)間在各小區(qū)間 上對(duì)其兩邊進(jìn)行積分，即上對(duì)其兩邊進(jìn)行積分，即 yxfy, 1, kkxx 111, 1 , 0,kkkkxxxxnkdxxyxfdxy 011,yaydxxyxfxyxykkxxkk如用矩形數(shù)值積分公式可得：如用矩形數(shù)值積分公式可得： 1, 1 , 0,1

7、0nkyxhfyyyaykkkk以上三種方法推導(dǎo)出同一個(gè)數(shù)值求解公式以上三種方法推導(dǎo)出同一個(gè)數(shù)值求解公式: :這個(gè)數(shù)值公式稱為這個(gè)數(shù)值公式稱為歐拉歐拉(Euler)(Euler)公式。公式。 1, 1 , 0,10nkyxhfyyyaykkkk6.1 歐拉方法歐拉方法一、一、 歐拉格式：歐拉格式：x0 x1向前差商近似導(dǎo)數(shù)向前差商近似導(dǎo)數(shù)hxyxyxy)()()(010 ),()()()(000001yxfhyxyhxyxy 1y記為記為歐拉公式幾何意義歐拉公式幾何意義 用一條通過初始點(diǎn)的折用一條通過初始點(diǎn)的折線近似表示解曲線線近似表示解曲線 ,亦稱為亦稱為歐歐拉折線法拉折線法 ，或稱為或稱為

8、矩形法。矩形法。)1,., 0(),(1 nkyxfhyykkkk一般形式一般形式1 1、顯式歐拉公式、顯式歐拉公式在假設(shè)在假設(shè) yk = y(xk)，即第，即第 k 步計(jì)算是精確的前提下，考步計(jì)算是精確的前提下，考慮的截?cái)嗾`差慮的截?cái)嗾`差 Rk = y(xk+1) yk+1 稱為稱為局部截?cái)嗾`差局部截?cái)嗾`差 。定義定義 若某算法的局部截?cái)嗾`差為若某算法的局部截?cái)嗾`差為O(hp+1)，則稱該算法有，則稱該算法有p 階精度。階精度。定義定義 歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差：),()()()()(2112kkyxhfyyxyhxyyxyRkhkkkki )(22 yh 歐拉法具有歐拉法

9、具有 1 階精度。階精度。局部截?cái)嗾`差和階數(shù)局部截?cái)嗾`差和階數(shù)1 kkxx 2、隱式歐拉格式、隱式歐拉格式向后差商近似導(dǎo)數(shù)向后差商近似導(dǎo)數(shù)hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy )1,., 0(),(111 nkyxfhyykiik由于未知數(shù)由于未知數(shù) yk+1 同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故稱為稱為隱式隱式 歐拉公式，而前者稱為歐拉公式，而前者稱為顯式顯式 歐拉公式。歐拉公式。一般先用顯式計(jì)算一個(gè)初值，再一般先用顯式計(jì)算一個(gè)初值，再迭代迭代求解。求解。 隱式隱式歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差：

10、11)( kkkyxyR)(22 yh 即隱式歐拉公式具有即隱式歐拉公式具有 1 階精度。階精度。二、兩步歐拉格式（中點(diǎn)公式）二、兩步歐拉格式（中點(diǎn)公式）中心差商近似導(dǎo)數(shù)中心差商近似導(dǎo)數(shù)hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 nkyxfhyykkkk假設(shè)假設(shè) ,則可以導(dǎo)出則可以導(dǎo)出即兩步歐拉格式具有即兩步歐拉格式具有 2 階精度。階精度。)(),(11kkkkxyyxyy 311()()kkkRy xyO h 該方法需要該方法需要2個(gè)初值個(gè)初值 y0和和 y1來啟動(dòng)遞推過程，這樣的算法稱來啟動(dòng)遞推過程，這樣的算法

11、稱為為雙步法雙步法。三、三、 梯形公式梯形公式 顯、隱式兩種算法的顯、隱式兩種算法的平均平均)1,., 0(),(),(2111 nkyxfyxfhyykkkkkk注：注：有局部截?cái)嗾`差有局部截?cái)嗾`差 ， 即梯形公式具有即梯形公式具有2 階精度，比歐拉方法有了進(jìn)步。階精度，比歐拉方法有了進(jìn)步。但注意到該公式是但注意到該公式是隱式隱式公式，計(jì)算時(shí)不得不用到公式，計(jì)算時(shí)不得不用到迭代法迭代法，不易求解。，不易求解。)()(311hOyxyRkkk 對(duì)歐拉法進(jìn)行改進(jìn)，用梯形公式計(jì)算右側(cè)積分，即對(duì)歐拉法進(jìn)行改進(jìn)，用梯形公式計(jì)算右側(cè)積分，即 11,2,1 kkkkxxyxfyxfhdxyxfkk 2 ,

12、 1 , 0,),(),(2),()0(11)1(1)0(1 ikyxfyxfhyyyxhfyykkkkkikkkkk計(jì)算計(jì)算公式公式梯形格式算法計(jì)算步驟：梯形格式算法計(jì)算步驟： 先用先用(1)式計(jì)算出式計(jì)算出 處處 。1kx (0)1ky 再用再用(2)式反復(fù)進(jìn)行迭代，得到式反復(fù)進(jìn)行迭代，得到(1)(2)11,kkyy(0)1(1)(0)111(,)(,)(,)2kkkkikkkkkkyyhf xyhyyf xyf xy 計(jì)算計(jì)算公式公式-(1)-(2)類似地得到類似地得到 用用 控制迭代次數(shù)，控制迭代次數(shù)， 為允許誤差。為允許誤差。把滿足誤差要求的把滿足誤差要求的 作為作為 的近似值。的近

13、似值。 (1)( )11iikkyy (1)1iky 1ky x 23,kkyy方方 法法 顯式歐拉法顯式歐拉法 隱式歐拉法隱式歐拉法 梯形公式梯形公式 中點(diǎn)公式中點(diǎn)公式 簡(jiǎn)單簡(jiǎn)單 精度低精度低 穩(wěn)定性最好穩(wěn)定性最好 精度低精度低, 計(jì)算量大計(jì)算量大 精度提高精度提高 計(jì)算量大計(jì)算量大 精度提高精度提高, 顯式顯式 多一個(gè)初值多一個(gè)初值, 可能影響精度可能影響精度 不同方法比較不同方法比較四、四、改進(jìn)歐拉法（預(yù)報(bào)改進(jìn)歐拉法（預(yù)報(bào)-校正法）校正法）Step 1: 先用先用顯式顯式歐拉公式作歐拉公式作預(yù)報(bào)預(yù)報(bào)，算出，算出),(1iiiiyxfhyy Step 2: 再將再將 代入代入隱式隱式梯形公

14、式的右邊作梯形公式的右邊作校正校正，得到，得到1 iy),(),(2111 iiiiiiyxfyxfhyy )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii它可表示為嵌套形式它可表示為嵌套形式表示為平均化形式表示為平均化形式 pcipiiciiipyyyyxhfyyyxhfyy21,11此法稱為此法稱為預(yù)報(bào)預(yù)報(bào)-校正法，校正法，是是顯式算法顯式算法。注：注：可以證明該算法具有可以證明該算法具有 2 階精度，同時(shí)可以看到它是個(gè)階精度，同時(shí)可以看到它是個(gè)單單步步遞推格式（只迭代一次）遞推格式（只迭代一次） ，比隱式梯形公式的迭代，比隱式梯形公式的迭代求解過程求解

15、過程簡(jiǎn)單簡(jiǎn)單。它的。它的穩(wěn)定性高穩(wěn)定性高于顯式歐拉法。于顯式歐拉法。腳標(biāo)用腳標(biāo)用 i舉例：進(jìn)行比較。進(jìn)行比較。精確解精確解并與并與法求解法求解法和改進(jìn)法和改進(jìn)試分別用試分別用設(shè)初值問題設(shè)初值問題例例xyyyxydxdy21,EulerEuler1)0(2 計(jì)算結(jié)果如下表所示：計(jì)算結(jié)果如下表所示：法：法：改進(jìn)的改進(jìn)的法：法：上結(jié)果，此時(shí)上結(jié)果，此時(shí)計(jì)算計(jì)算解：取解：取 ,.)2 , 1 , 0()(21)2()2(Euler,.)2 , 1 , 0()2( 1 . 01 , 0, 1 . 0111iyyyyxyhyyyxyhyyiyxyyyEulerxhcpipipiciiiipiiiiixEu

16、ler法y改進(jìn)的Euler法y精確解01.0000001.0000001.0000000.11.0000001.0959091.0954450.21.1918181.1840971.1832160.31.2774381.2662011.2649110.41.3582131.3433601.3416410.51.4351331.4164021.4142140.61.5089661.4859561.4832400.71.5803381.5525141.5491930.81.6497831.6164751.6124520.91.7177791.6781661.6733201.01.7847701.7

17、378671.7320516.2 龍格龍格 - 庫塔法庫塔法一、一、 泰勒級(jí)數(shù)法泰勒級(jí)數(shù)法 龍格龍格庫塔庫塔(Runge-Kut ta)法法(簡(jiǎn)稱為簡(jiǎn)稱為R-K方法方法)是一類高精是一類高精度的一步法，這類方法與泰勒級(jí)數(shù)法有著密度的一步法，這類方法與泰勒級(jí)數(shù)法有著密 切的關(guān)系。切的關(guān)系。 設(shè)有初值問題設(shè)有初值問題 00)()(,()(yxyxyxfxy由由 泰勒展開式泰勒展開式 1)(21! 2)( kkkkkkkkhxykhxyhxhyxyxy 從理論上講，只要解從理論上講，只要解y(x)有任意階導(dǎo)數(shù)，泰勒展開方有任意階導(dǎo)數(shù)，泰勒展開方法就可以構(gòu)造法就可以構(gòu)造任意階任意階求求yk+1公式。但

18、由于計(jì)算這些導(dǎo)數(shù)是公式。但由于計(jì)算這些導(dǎo)數(shù)是非常復(fù)雜的，所以這種方法實(shí)際上不能用來解初值問題。非常復(fù)雜的，所以這種方法實(shí)際上不能用來解初值問題。 設(shè)有初值問題設(shè)有初值問題 00)()(,()(yxyxyxfxy二、二、龍格庫塔法的基本思路龍格庫塔法的基本思路)(),()(hxyhxhfxykkk 1)(,()()(1kkxxkkdxxyxfxyxy等價(jià)于：等價(jià)于：（積分中值定理）（積分中值定理） R-K方法基本思想：方法基本思想：用用 在幾個(gè)不同點(diǎn)的加權(quán)平均值在幾個(gè)不同點(diǎn)的加權(quán)平均值（線性組合）來代替準(zhǔn)確的（線性組合）來代替準(zhǔn)確的 的值，構(gòu)的值，構(gòu)造近似公式。再把近似公式與解的泰勒展開造近似公

19、式。再把近似公式與解的泰勒展開 式進(jìn)行比較，式進(jìn)行比較，使前面的若干項(xiàng)相同，從而使近似公式達(dá)到一定的階數(shù)。使前面的若干項(xiàng)相同，從而使近似公式達(dá)到一定的階數(shù)。 這樣龍格庫塔法保留了泰勒級(jí)數(shù)展開法的高階局部截這樣龍格庫塔法保留了泰勒級(jí)數(shù)展開法的高階局部截?cái)嗾`差，又避免了高階導(dǎo)數(shù)的計(jì)算。斷誤差，又避免了高階導(dǎo)數(shù)的計(jì)算。 我們先分析歐拉法我們先分析歐拉法 與預(yù)估與預(yù)估校正法。校正法。),(yxf)(,(hxyhxfkk 112111121(,)(,)(,)(,)1()()22pkkkkkckkpkkkkkpcyyhf xyKf xyyyhf xyKf xyhKhyyKKyyy 改改進(jìn)進(jìn)的的歐歐拉拉式式

20、在在公公中中1111131(,)(,)(,)(,),()kkkkkkkkkkkkxyyf xyxyhKyf xyhKxxO h 則則是是用用點(diǎn)點(diǎn)處處的的斜斜率率和和由由此此點(diǎn)點(diǎn)處處信信息息預(yù)預(yù)估估的的點(diǎn)點(diǎn)處處的的斜斜率率的的算算術(shù)術(shù)平平均均值值來來近近似似代代替替區(qū)區(qū)間間上上的的平平均均斜斜率率，局局部部截截?cái)鄶嗾`誤差差為為。0121( )(,)(,)(,),()kkkkkkkkkkky ayyyhf xyxyyf xyxxO h 在在歐歐拉拉公公式式中中僅僅用用一一個(gè)個(gè)點(diǎn)點(diǎn)處處的的斜斜率率來來近近似似代代替替區(qū)區(qū)間間上上的的平平均均斜斜率率，局局部部截截?cái)鄶嗾`誤差差為為。所所以以如如果果在在

21、區(qū)區(qū)間間上上多多預(yù)預(yù)估估幾幾個(gè)個(gè)點(diǎn)點(diǎn)的的斜斜率率值值，再再將將它它們們的的線線性性組組合合作作為為平平均均斜斜率率的的近近似似值值，則則就就有有可可能能構(gòu)構(gòu)造造出出精精度度更更高高的的計(jì)計(jì)算算格格式式。推廣推廣1(, )kkkkyyhxyh ( , )f x y 其其中中 是是用用在在一一些些點(diǎn)點(diǎn)上上值值的的線線性性組組合合來來構(gòu)構(gòu)成成這種單步法稱為這種單步法稱為Runge-KuttaRunge-Kutta方法方法, ,簡(jiǎn)記為簡(jiǎn)記為R-KR-K公式公式. .1(, )rkkiiixyhc k K Ki i為某些點(diǎn)上的斜率，或?yàn)槟承c(diǎn)上的斜率，或f(x,yf(x,y) )在某些點(diǎn)上的值。在某些點(diǎn)

22、上的值。,Runge-Kutta.若若 是是由由 個(gè)個(gè) 值值線線性性組組合合構(gòu)構(gòu)成成 則則稱稱線線性性階階方方法法RRf 三、三、 二階龍格二階龍格 - 庫塔法庫塔法目標(biāo)目標(biāo)：建立高精度的單步遞推格式。建立高精度的單步遞推格式。單步遞推法的單步遞推法的基本思想基本思想是從是從 ( xi , yi ) 點(diǎn)出發(fā)，以點(diǎn)出發(fā)，以某一斜某一斜率率沿直線達(dá)到沿直線達(dá)到 ( xi+1 , yi+1 ) 點(diǎn)。歐拉法及其各種變形所點(diǎn)。歐拉法及其各種變形所能達(dá)到的最高精度為能達(dá)到的最高精度為2階階。 考察改進(jìn)的歐拉法，可以將其改寫為：考察改進(jìn)的歐拉法，可以將其改寫為：),(),(2121121211hKyhxfK

23、yxfKKKhyyiiiiii 斜率斜率一定取一定取K1 K2 的的平均值平均值嗎？嗎？ 步長(zhǎng)一定是一個(gè)步長(zhǎng)一定是一個(gè)h 嗎？嗎？ 腳標(biāo)用腳標(biāo)用 i首先希望能確定系數(shù)首先希望能確定系數(shù) 1、 2、p，使得到的算法格式有，使得到的算法格式有2階階精度，即在精度，即在 的前提假設(shè)下，使得的前提假設(shè)下，使得 )(iixyy )()(311hOyxyRiii Step 1: 將將 K2 在在 ( xi , yi ) 點(diǎn)作點(diǎn)作 Taylor 展開展開)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 將改進(jìn)歐拉法

24、推廣為：將改進(jìn)歐拉法推廣為：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2: 將將 K2 代入代入yi+1表達(dá)式，得到表達(dá)式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii Step 3: 將將 yi+1 與與 y( xi+1 ) 在在 xi 點(diǎn)的點(diǎn)的泰勒泰勒展開作比較展開作比較)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhx

25、yhxyxyiiii 要求要求 ，則必須有：，則必須有：)()(311hOyxyRiii21,1221 p 這里有這里有 個(gè)未知個(gè)未知數(shù)，數(shù)， 個(gè)方程。個(gè)方程。32存在存在無窮多個(gè)解無窮多個(gè)解。所有滿足上式的格式統(tǒng)稱為。所有滿足上式的格式統(tǒng)稱為2階龍格階龍格 - 庫庫塔格式塔格式。21, 121 p注意到，注意到， 就是改進(jìn)的歐拉法。就是改進(jìn)的歐拉法。 Q: 為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均為待定系數(shù)

26、，確定這些系數(shù)的步驟與前面相似。其均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似。其解不唯一。解不唯一。).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用為四階最常用為四階經(jīng)典龍格經(jīng)典龍格-庫塔法庫塔法四階經(jīng)典龍格四階經(jīng)典龍格-庫塔法公式庫塔法公式),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 四、四、 四階龍格四階龍格 - 庫塔法庫塔法用四個(gè)用

27、四個(gè)f f函數(shù)值的線性組合得到四階函數(shù)值的線性組合得到四階龍格龍格 - - 庫塔法庫塔法。經(jīng)典龍格經(jīng)典龍格-庫塔法公式具有四階精度，因此可取大步長(zhǎng)。庫塔法公式具有四階精度，因此可取大步長(zhǎng)。注：注： 龍格龍格-庫塔法庫塔法的主要運(yùn)算在于計(jì)算的主要運(yùn)算在于計(jì)算 Ki 的值，即計(jì)算的值，即計(jì)算 f 的的值。值。Butcher 于于1965年給出了計(jì)算量與可達(dá)到的最高精年給出了計(jì)算量與可達(dá)到的最高精度階數(shù)的關(guān)系：度階數(shù)的關(guān)系：753可達(dá)到的最高精度可達(dá)到的最高精度642每步須算每步須算Ki 的個(gè)數(shù)的個(gè)數(shù)2345642n 8n高于四階時(shí)每步計(jì)算量增加較多，但精度提高不快，高于四階時(shí)每步計(jì)算量增加較多，但

28、精度提高不快，因此使用的比較少。因此使用的比較少。由于龍格由于龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對(duì)于光滑性不太好的解，最好解函數(shù)的光滑性影響。對(duì)于光滑性不太好的解，最好采用采用低階算法低階算法而將步長(zhǎng)而將步長(zhǎng)h 取小取小。Lab 15. Runge-Kutta (Order Four) Use Runge-Kutta method of order four to approximate the solution of the initial-value problem, , and (1)You are supposed t

29、o write a functionvoid RK4 ( double (*f)( ), double a, double b, double y0, int n, FILE *outfile)to approximate the solution of Problem (1) with y= f, x in a, b, and the initial value of y being y0. Output the approximating values of y on the n+1 equally spaced grid points from a to b to outfile.Inp

30、utThere is no input file. Instead, you must hand in your function in a *.h file. The rule of naming the *.h file is the same as that of naming the *.c or *.cpp files.Output ( represents a space) For each test case, you are supposed to print n+1 lines, and each line is in the following format: fprint

31、f( outfile, “%8.4f%12.8en”, x, y );),(yxfy ,bax 0)(yay Sample Judge ProgramSample Judge Program#include #include #include 98115001_15.h double f0(double x, double y) return (y x*x+1.0); void main( ) FILE *outfile = fopen(out.txt, w);int n;double a, b, y0; a = 0.0; b = 2.0; y0 = 0.50; n = 10;RK4(f0,

32、a, b, y0, n, outfile);fprintf(outfile, n);fclose(outfile);Sample Output Sample Output ( ( represents a space) represents a space) 0.00005.00000000e 0010.20008.29293333e 0010.40001.21407621e+0000.60001.64892202e+0000.80002.12720268e+0001.00002.64082269e+0001.20003.17989417e+0001.40003.73234007e+0001.

33、60004.28340950e+0001.80004.81508569e+0002.00005.30536300e+0003 收斂性與穩(wěn)定性收斂性與穩(wěn)定性 /* Convergency and Stability */ 收斂性收斂性 /* Convergency */定義定義 若某算法對(duì)于任意固定的若某算法對(duì)于任意固定的 x = xi = x0 + i h，當(dāng)，當(dāng) h0 ( 同時(shí)同時(shí) i ) 時(shí)有時(shí)有 yi y( xi )，則稱該算法是，則稱該算法是收斂收斂的。的。 例：例：就初值問題就初值問題 考察歐拉顯式格式的收斂性?？疾鞖W拉顯式格式的收斂性。 0)0(yyyy 解：解：該問題的精確解為該

34、問題的精確解為 xeyxy 0)( 歐拉公式為歐拉公式為iiiiyhyhyy)1 (1 0)1 (yhyii 對(duì)任意固定的對(duì)任意固定的 x = xi = i h ，有，有iixhhxihyhyy )1()1(/10/0 ehhh /10)1(lim)(0ixxyeyi 穩(wěn)定性穩(wěn)定性 /* Stability */例：例：考察初值問題考察初值問題 在區(qū)間在區(qū)間0, 0.5上的解。上的解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計(jì)算數(shù)值解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計(jì)算數(shù)值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精確解精確解改進(jìn)歐拉法改進(jìn)歐拉法 歐歐拉隱式拉隱式歐拉顯式歐拉顯式 節(jié)點(diǎn)節(jié)點(diǎn) xixey30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 10 21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4