計(jì)算方法Chapter03-常微分方程的數(shù)值解法

1、三常微分方程數(shù)值解法2常微分方程數(shù)值解法引言（常微分方程數(shù)值解法概述）引言（常微分方程數(shù)值解法概述）顯式歐拉法、隱式歐拉法、二步歐拉法顯式歐拉法、隱式歐拉法、二步歐拉法局部截?cái)嗾`差與精度局部截?cái)嗾`差與精度改進(jìn)的歐拉方法改進(jìn)的歐拉方法龍格龍格-庫塔方法庫塔方法收斂性與穩(wěn)定性簡(jiǎn)述收斂性與穩(wěn)定性簡(jiǎn)述一階常微分方程組與高階常微分方程一階常微分方程組與高階常微分方程3引言引言一階常微分方程初值問題：一階常微分方程初值問題：00( , )()yf x yy xy 微分方程微分方程初始條件初始條件定理：定理：若若 f (x, y) 在某閉區(qū)域在某閉區(qū)域 R ： 00|,|(0,0)xxayybab 上連續(xù)，

2、且在上連續(xù)，且在 R 域內(nèi)滿足李普希茲域內(nèi)滿足李普希茲 (Lipschitz) 條件，條件，即存在正數(shù)即存在正數(shù) L，使得對(duì)于，使得對(duì)于 R 域內(nèi)的任意兩值域內(nèi)的任意兩值 y1, y2，下，下列不等式成立：列不等式成立： 1221|( ,)( ,)|f x yf x yL yy 則上述初值問題的連續(xù)可微的解則上述初值問題的連續(xù)可微的解 y(x) 存在并且唯一。存在并且唯一。4引言（續(xù)）引言（續(xù)）實(shí)際生產(chǎn)與科研中，除少數(shù)簡(jiǎn)單情況能獲得初值問題實(shí)際生產(chǎn)與科研中，除少數(shù)簡(jiǎn)單情況能獲得初值問題的初等解（用初等函數(shù)表示的解）外，絕大多數(shù)情況的初等解（用初等函數(shù)表示的解）外，絕大多數(shù)情況下是求不出初等解的

3、。下是求不出初等解的。有些初值問題即便有初等解，也往往由于形式過于復(fù)有些初值問題即便有初等解，也往往由于形式過于復(fù)雜而不便處理。雜而不便處理。實(shí)用的方法是在計(jì)算機(jī)上進(jìn)行數(shù)值求解：即不直接求實(shí)用的方法是在計(jì)算機(jī)上進(jìn)行數(shù)值求解：即不直接求 y(x) 的顯式解，而是在解所存在的區(qū)間上，求得一系的顯式解，而是在解所存在的區(qū)間上，求得一系列點(diǎn)列點(diǎn) xn (n 0, 1, 2, ) 上解的近似值。上解的近似值。5歐拉歐拉（Euler）方法方法方法一方法一化導(dǎo)數(shù)為差商的方法化導(dǎo)數(shù)為差商的方法10()()()()()limnnnnnhy xhy xy xy xy xhh 由于在逐步求解的過程中，由于在逐步求解

4、的過程中，y(xn) 的準(zhǔn)確值無法求解的準(zhǔn)確值無法求解出來，因此用其近似值代替。出來，因此用其近似值代替。為避免混淆，以下學(xué)習(xí)簡(jiǎn)記：為避免混淆，以下學(xué)習(xí)簡(jiǎn)記：y(xn)：待求函數(shù)：待求函數(shù) y(x) 在在 xn 處的精確函數(shù)值處的精確函數(shù)值yn ：待求函數(shù)：待求函數(shù) y(x) 在在 xn 處的近似函數(shù)值處的近似函數(shù)值6代入初值問題表達(dá)式可得：代入初值問題表達(dá)式可得：根據(jù)根據(jù) y0 可以一步步計(jì)算出函數(shù)可以一步步計(jì)算出函數(shù) y y(x) 在在 x1, x2, x3 x4, 上的近似值上的近似值 y1, y2, y3, y4 , 常微分方程數(shù)值解是一組離散的函數(shù)值數(shù)據(jù)，它的常微分方程數(shù)值解是一組離

5、散的函數(shù)值數(shù)據(jù)，它的精確表達(dá)式很難求解得到，但可以進(jìn)行插值計(jì)算后精確表達(dá)式很難求解得到，但可以進(jìn)行插值計(jì)算后用插值函數(shù)逼近用插值函數(shù)逼近 y(x)10()()()()()limnnnnnhy xhy xy xy xy xhh 100(,)0,1,2,()nnnnyyhf xynyy x 1(,)nnnnyyf xyh 7歐拉方法（續(xù)）歐拉方法（續(xù)）方法二方法二數(shù)值積分法數(shù)值積分法1111d()() , ( )d() , (), ()nnnnxxnnxxnnnnnnyxy xy xf x y xxxxf xy xhf xy x 同樣以近似值同樣以近似值 yn 代替精確值代替精確值 y(xn) 可

6、得可得：100(,)()nnnnyyhf xyyy x 將微分方程將微分方程 y f (x, y) 在區(qū)間在區(qū)間 xn, xn+1 上積分：上積分：8歐拉方法的幾何意義歐拉方法的幾何意義xy00 x1x2x3x4x5x6x0P1P2P3P4P5P6P( )yy x9隱式歐拉法隱式歐拉法在數(shù)值積分法推導(dǎo)中，積分的近似值取為積分區(qū)間寬在數(shù)值積分法推導(dǎo)中，積分的近似值取為積分區(qū)間寬度與右端點(diǎn)處的函數(shù)值乘積，即：度與右端點(diǎn)處的函數(shù)值乘積，即：這樣便得到了隱式歐拉法：這樣便得到了隱式歐拉法：11100(,)()nnnnyyh f xyyy x 含有未知含有未知的函數(shù)值的函數(shù)值隱式歐拉法沒有顯式歐拉法方便

7、隱式歐拉法沒有顯式歐拉法方便11111111d()() , ( )d() , (), ()nnnnxxnnxxnnnnnnyxy xy xf x y xxxxf xy xhf xy x 10二步歐拉法二步歐拉法在數(shù)值積分法推導(dǎo)中，積分區(qū)間寬度選為兩步步長，在數(shù)值積分法推導(dǎo)中，積分區(qū)間寬度選為兩步步長，即積分區(qū)間為：即積分區(qū)間為：xn 1, xn 1，則：，則：11111111d()() , ( )d() , ()2, ()nnnnxxnnxxnnnnnnyxy xy xf x y xxxxf xy xh f xy x 以以 y(x) 在在 xn 1, xn 上的近似值上的近似值代替精確值代替精

8、確值可得可得：11002(,)()nnnnyyhf xyyy x 需要前兩步需要前兩步的計(jì)算結(jié)果的計(jì)算結(jié)果中矩形公式中矩形公式11梯形公式歐拉法梯形公式歐拉法在數(shù)值積分法中，如果用梯形公式近似計(jì)算在數(shù)值積分法中，如果用梯形公式近似計(jì)算 f (x, y) 在區(qū)間在區(qū)間 xn, xn+1 上的積分，即：上的積分，即：用近似值代替精確值可得梯形公式歐拉法：用近似值代替精確值可得梯形公式歐拉法：111 (,)(,)2nnnnnnhyyf xyf xy 上式右端出現(xiàn)了未知項(xiàng)，可見梯形法是隱式歐拉法上式右端出現(xiàn)了未知項(xiàng)，可見梯形法是隱式歐拉法的一種；實(shí)際上，梯形公式歐拉法是顯式歐拉法與的一種；實(shí)際上，梯形

9、公式歐拉法是顯式歐拉法與隱式歐拉法的隱式歐拉法的算術(shù)平均算術(shù)平均。 111111, (), () , ( )d()2, (), ()2nnxnnnnnnxnnnnf xy xf xy xf x y xxxxhf xy xf xy x 12例例用顯式歐拉法、隱式歐拉法、梯形法求解初值問題：用顯式歐拉法、隱式歐拉法、梯形法求解初值問題：1(0)1yyxy 取取 h 0.1，計(jì)算到，計(jì)算到 x 0.5，并與精確解進(jìn)行比較，并與精確解進(jìn)行比較解：由已知條件可得：解：由已知條件可得：h 0.1，x0 0， y0 1， f (x, y) y x 1顯式歐拉法：顯式歐拉法：1(1)(1)0.90.10.1n

10、nnnnnnnyyhyxh yhxhyx 13例：（續(xù)）例：（續(xù)）隱式歐拉法：隱式歐拉法：111(1)nnnnyyhyx 化簡(jiǎn)得：化簡(jiǎn)得：11(1)(1)(1)nnnnnh yyh xyh xh 11(1)(0.10.11) 1.1(1)nnnnnyyh xhyxh 111(1)(1)2nnnnnnhyyyxyx 梯形公式歐拉法：梯形公式歐拉法：111(2)1(2)(2)(2)nnnnhyhyhxx 11(2)(2) (2)(1.90.20.21) 2.1nnnnnnyh yh xxhyx 14計(jì)算結(jié)果：計(jì)算結(jié)果：xn顯式法顯式法 yn隱式法隱式法 yn梯形法梯形法 yn精確解精確解 y (x

11、n)0.011110.11.0000001.0090911.0047621.0048370.21.0100001.0264461.0185941.0197310.31.0290001.0513151.0406331.0408180.41.0561001.0830141.0700971.0703200.51.0904901.1209221.1062781.106531( )exy xx 本題的精確解為：本題的精確解為：10.90.10.1nnnyyx 1(0.10.11) 1.1nnnyyx 1(1.90.20.21) 2.1nnnyyx 顯式法顯式法隱式法隱式法梯形法梯形法15局部截?cái)嗾`差局部

12、截?cái)嗾`差為了簡(jiǎn)化分析某常微分方程數(shù)值算法的誤差，現(xiàn)假為了簡(jiǎn)化分析某常微分方程數(shù)值算法的誤差，現(xiàn)假設(shè)設(shè) yn y(xn)，即，即在前一步在前一步 yn 準(zhǔn)確的前提下準(zhǔn)確的前提下，估計(jì)：，估計(jì)：111()nnnTy xy 稱上述誤差稱上述誤差 Tn 1 為該常微分方程數(shù)值算法的為該常微分方程數(shù)值算法的局部截局部截?cái)嗾`差斷誤差如果某個(gè)常微分方程數(shù)值算法的局部截?cái)嗾`差可表如果某個(gè)常微分方程數(shù)值算法的局部截?cái)嗾`差可表示為示為 O(h p 1)，則稱該數(shù)值算法的精度是，則稱該數(shù)值算法的精度是 p 階階歐拉法的精度為一階；二步歐拉法的精度為二階；歐拉法的精度為一階；二步歐拉法的精度為二階；梯形公式歐拉法的精

13、度為二階。梯形公式歐拉法的精度為二階。16泰勒展開法泰勒展開法如果初值問題中的如果初值問題中的 f (x, y) 充分可微，則可將充分可微，則可將 y(xn 1) 在點(diǎn)在點(diǎn) xn 處展開：處展開：21()()()()()2!nnnnnhy xy xhy xhy xyx 如果只保留線性項(xiàng)，忽略如果只保留線性項(xiàng)，忽略 h2 及以上各項(xiàng)，則：及以上各項(xiàng)，則：11()()()nnnny xy xhy xy (), ()nnny xf xy x ny1(,)nnnnyyhf xy 顯式歐拉公式顯式歐拉公式17局部截?cái)嗾`差的分析局部截?cái)嗾`差的分析利用泰勒公式展開，比較各算法與展開式的前幾項(xiàng)利用泰勒公式展開

14、，比較各算法與展開式的前幾項(xiàng)將將 y(xn 1) 在在 xn 點(diǎn)處用泰勒公式展開：點(diǎn)處用泰勒公式展開：211( )()()()()2!(,)nnnnnnyy xy xhy xhy xhxx 22111( )()()2!nnnyTy xyhO h 顯式歐拉法的局部截?cái)嗾`差：顯式歐拉法的局部截?cái)嗾`差：1(), ()()()nnnnnnyy xhf xy xy xhy x 1(,)nnnnyyhf xy 歐拉法歐拉法()nnyy x 令令：1 階精度階精度18補(bǔ)充：二元函數(shù)微分中值定理補(bǔ)充：二元函數(shù)微分中值定理1100(,)(,)f xyf xy 1001010010()()()()xyxxfxxx

15、yyfyyy 01 0000(,)(,)f xh ykf xy 00()()xyhfxhkfyk 000000(,)(,)()()xyf xh ykf xyhfxhkfyk 0000(,)(,)f xyhkf xh ykxy 191111111()(,),(, )()nnnynnnnyf xyf xxxfxyy 111()()( )(, )nnnynnyy xh y xO hfxT 21()()()()nnny xy xhy xO h y(xn 1) 在在 xn 點(diǎn)處展開：點(diǎn)處展開：211111()(, )()nnnynnTy xyhfxTO h 111(,)nnnnyyhf xy 隱式歐拉法

16、：隱式歐拉法：1 階精度階精度()ny x2111()1(, )nynTO hhfx 11 (),nny xy 111()(, )nynny xfxT 11()( )(, )nynny xO hfxT 211()()(, )()nnynny xhy xhfxTO h 20111()2, ()()2()nnnnnnyy xhf xy xy xhy x 分別將分別將 y(xn 1), y(xn 1) 在在 xn 點(diǎn)處用泰勒公式展開：點(diǎn)處用泰勒公式展開：231231()()()()()2!()()()()()()2!nnnnnnnnyxy xy xhy xhO hyxy xy xhy xhO h 1

17、11()nnnTy xy 二步歐拉法的局部截?cái)嗾`差：二步歐拉法的局部截?cái)嗾`差：112(,)nnnnyyhf xy 二步歐拉法：二步歐拉法：2 階精度階精度1()ny x ()ny x311()()2()()nnny xy xhy xO h 211111111111(,),(, )()()()(, )nnnynnnnynnnf xyf xfxyy xy xfTy xx 111 (,)(,)2nnnnnnhyyf xyf xy 梯形公式歐拉法：梯形公式歐拉法：()ny x2111()2()()()(, )2nnnnynnhyy xy xhyxO hfxT y(xn 1) 在在 xn 點(diǎn)處展開：點(diǎn)處

18、展開：231()()()()()2!nnnnhy xy xhy xyxO h 3111(, )()2nynnhTfxTO h 31211()1(, )nhynTO hfx 2 階精度階精度11 (),nny xy 2311()()()(, )()22nnnynnhhy xhy xyxfxTO h 211()()()(, )nnynny xhyxO hfxT 22各種歐拉法的比較各種歐拉法的比較方法方法精度精度評(píng)述評(píng)述顯式歐拉法顯式歐拉法1最簡(jiǎn)單，精度低最簡(jiǎn)單，精度低隱式歐拉法隱式歐拉法1不便計(jì)算，穩(wěn)定性好不便計(jì)算，穩(wěn)定性好二步歐拉法二步歐拉法2需要兩步初值，且第需要兩步初值，且第 2 個(gè)初值只

19、能由其個(gè)初值只能由其它方法給出，可能對(duì)后面的遞推精度有它方法給出，可能對(duì)后面的遞推精度有影響影響梯形公式梯形公式歐歐 拉拉 法法2精度有所提高，但為隱式，需要迭代求精度有所提高，但為隱式，需要迭代求解，計(jì)算量大解，計(jì)算量大23改進(jìn)的歐拉法改進(jìn)的歐拉法從上述例子可以看到，梯形法由于具有二階精度，其從上述例子可以看到，梯形法由于具有二階精度，其局部截?cái)嗾`差比顯式歐拉法和隱式歐拉法小，但梯形局部截?cái)嗾`差比顯式歐拉法和隱式歐拉法小，但梯形法實(shí)質(zhì)上是一種隱式算法法實(shí)質(zhì)上是一種隱式算法顯式歐拉法是一個(gè)顯式算法，雖然計(jì)算量較小，但是顯式歐拉法是一個(gè)顯式算法，雖然計(jì)算量較小，但是精度不高精度不高綜合兩種方法的

20、長處，可以先用顯式歐拉法求出綜合兩種方法的長處，可以先用顯式歐拉法求出 y(xn+1) 的一個(gè)粗略近似值，然后用它代入梯形法公式的一個(gè)粗略近似值，然后用它代入梯形法公式的右端，用梯形法計(jì)算的右端，用梯形法計(jì)算 y(xn+1) 的較為精確的近似值。的較為精確的近似值。1ny 預(yù)預(yù)報(bào)報(bào)值值1ny 校校正正值值24改進(jìn)的歐拉法（續(xù)）改進(jìn)的歐拉法（續(xù)）按照上述思想，可以建立如下預(yù)報(bào)按照上述思想，可以建立如下預(yù)報(bào)-校正系統(tǒng)：校正系統(tǒng)：1(,)nnnnyyh f xy 預(yù)預(yù)報(bào)報(bào): :按以上兩式求解常微分方程的算法稱為改進(jìn)的歐拉按以上兩式求解常微分方程的算法稱為改進(jìn)的歐拉法，它還可以表示為：法，它還可以表示

21、為： 11(,),(,)2nnnnnnnnhyyf xyf xyh f xy 嵌套形式嵌套形式11(,)(,)() 2pnnncnnpnpcyyh f xyyyh f xyyyy 平均化形式平均化形式2 階精度階精度111 (,)(,)2nnnnnnhyyf xyf xy 校校正正: :25用改進(jìn)歐拉法求上例所述的初值問題并與歐拉法和用改進(jìn)歐拉法求上例所述的初值問題并與歐拉法和梯形法比較誤差的大小。梯形法比較誤差的大小。解解：采用改進(jìn)歐拉法的嵌套形式：采用改進(jìn)歐拉法的嵌套形式： 11(,),(,)2nnnnnnnnhyyf xyf xyh f xy 10.1,0.5(0)1yyxhxy 計(jì)計(jì)算

22、算到到 1(1)(,)12nnnnnnnhyyxyh f xyx (1)(1)12nnnnnnnhyyxyhyxxh ( 2)(2)122nnhhhhyxh 0.9050.0950.1nnyx 26計(jì)算結(jié)果計(jì)算結(jié)果xn改進(jìn)歐改進(jìn)歐拉法拉法 yn精確解精確解 y (xn) 誤誤 差差改進(jìn)歐拉法改進(jìn)歐拉法歐拉法歐拉法梯形法梯形法0.1 1.005000 1.0048371.6 10-44.8 10-37.5 10-50.2 1.019205 1.0197312.9 10-48.7 10-31.4 10-40.3 1.041218 1.0408184.0 10-41.2 10-21.9 10-40.

23、4 1.070802 1.0703204.8 10-41.4 10-22.2 10-40.5 1.107076 1.1065315.5 10-41.6 10-22.5 10-4可見，改進(jìn)歐拉法的誤差數(shù)量級(jí)與梯形法大致相同，可見，改進(jìn)歐拉法的誤差數(shù)量級(jí)與梯形法大致相同，而比歐拉法小得多。而比歐拉法小得多。 10.1,0.5(0)1yyxhxy 計(jì)計(jì)算算到到27改進(jìn)的歐拉法的意義改進(jìn)的歐拉法的意義11(,)(,)() 2pnnncnnpnpcyyh f xyyyh f xyyyy 改進(jìn)的歐拉法改進(jìn)的歐拉法的平均化形式的平均化形式1k1nyhk 2k12112222nnnhhkkyykkyh 11(

24、)()()( )(,)nnnnny xy xhy xhyxx y (xn 1) 在點(diǎn)在點(diǎn) xn 處的一階展開式為：處的一階展開式為：28改進(jìn)的歐拉法的幾何意義改進(jìn)的歐拉法的幾何意義2nyhkPnxny0 xyRh1211(,)(,)nnnnkf xykf xyhk 1nx1nyhkQ斜率: k21ny 1ny 1222nhhykkS斜率: k1291121211()/2( ,)(,)nnnnnnyyh kkkf x ykf xyhk 龍格龍格-庫塔庫塔(Runge-Kutta)方法方法11(,)(,)() 2pnnncnnpnpcyyh f xyyyh f xyyyy 改進(jìn)的歐拉法改進(jìn)的歐拉法

25、（2 2 階精度）階精度）11()()()( )() , ( )(,)nnnnnny xy xhy xhyy xh fyxx y (xn 1) 在點(diǎn)在點(diǎn) xn 處的一階泰勒展開式為：處的一階泰勒展開式為：*k1(,)nnnnyyhf xy 顯式歐拉法顯式歐拉法（1 1 階精度）階精度）111(,)nnnnyyhkkf xy 30龍格龍格-庫塔方法（續(xù)）庫塔方法（續(xù)）顯式歐拉法用一個(gè)點(diǎn)的值顯式歐拉法用一個(gè)點(diǎn)的值 k1 作為作為 k* 的近似值的近似值改進(jìn)的歐拉公式用二個(gè)點(diǎn)的值改進(jìn)的歐拉公式用二個(gè)點(diǎn)的值 k1 和和 k2 的平均值作為的平均值作為 k* 近似值；近似值；改進(jìn)的歐拉法比顯式歐拉法精度

26、高；改進(jìn)的歐拉法比顯式歐拉法精度高；在在 xn, xn 1 內(nèi)多預(yù)報(bào)幾個(gè)點(diǎn)的內(nèi)多預(yù)報(bào)幾個(gè)點(diǎn)的 ki 值，并用其加權(quán)平均值，并用其加權(quán)平均值作為值作為 k* 的近似值從而構(gòu)造出具有更高精度的計(jì)算公的近似值從而構(gòu)造出具有更高精度的計(jì)算公式，這就是式，這就是龍格龍格- -庫塔方法的基本思想。庫塔方法的基本思想。31二階龍格二階龍格-庫塔方法庫塔方法以以 k1 和和 k2 的加權(quán)平均來近似取代的加權(quán)平均來近似取代 k*12221121211()(,)(,)nnnnnncyyhkkkf xykf xh yhcakb 12221,ccab為為待待定定系系數(shù)數(shù)為分析局部截?cái)嗾`差，令為分析局部截?cái)嗾`差，令 y

27、n y(xn)，由泰勒公式得：，由泰勒公式得：231()()()()()()2!nnnnnhy xy xhy xhy xyxO h (), ()(,)nnnnny xf xy xf xy 123()()(,)(,)(,) (,)()2nnnnxnnynnnny xy xhf xyhfxyfxyf xyO h ()(,)(,)(,) (,)nnnxnnynnnnyxfxyfxyfxyf xy ny 32補(bǔ)充：二元泰勒展開式補(bǔ)充：二元泰勒展開式00000020000(,)(,)(,)1(,)2!1(,)!nf xh ykf xyhkf xyxyhkf xyxyhkf xynxy 20000200(

28、,)2(,)(,)xxxyyyh fxyhk fxyk fxy 000nniin inin i x xyyifC h kx y 0000(,)(,)xyhfxyk fxy 3322211(,)nnkf xa h yb hk 用二元泰勒公式展開用二元泰勒公式展開222211(,)(,)(,)()nnxnnynnkf xya hfxyb hk fxyO h 123()()(,)(,)(,) (,)()2nnnnxnnynnnny xy xhf xyhfxyfxyf xyO h 11223211(,)(,)(,)(,)()nnnnnnxnnynnyyh c f xycf xya h fxyb hk

29、fxyO h 將將 k1, k2 代入代入 中可得：中可得：11122()nnyyh c kc k 2122223221() (,)(,)(,) (,)()nnnxnnynnnnyh ccf xyh c a fxyh c b fxyf xyO h 122222213() (,)2(,)2(,) (,)2()nnnxnnynnnnyh ccf xyhc a fxyc b fxyf xyO h 1(,)nnkf xy 341122322221() (,)2(,)2(,) (,)()2nnnnxnnynnnnyyh ccf xyhc a fxyc b fxyf xyO h 123()()(,)(,)

30、(,) (,)()2nnnnxnnynnnny xy xhf xyhfxyfxyf xyO h 二階龍格二階龍格-庫塔方法（續(xù)）庫塔方法（續(xù)）122222112121ccc ac b 311()()nny xyO h 2 階精度階精度35四個(gè)未知變量，只有三四個(gè)未知變量，只有三個(gè)方程，有無窮多組解個(gè)方程，有無窮多組解每組解的構(gòu)成的龍格每組解的構(gòu)成的龍格-庫庫塔方法均為二階塔方法均為二階11122122211()(,)(,)nnnnnnyyh c kc kkf xykf xa h yb hk 122211 21ccab 取取二階龍格二階龍格- -庫塔方法庫塔方法即為改進(jìn)的歐拉方法即為改進(jìn)的歐拉方

31、法12221011 2ccab 取取12121(,),22nnnnnnyyhkkf xyhhkfxyk 變形的歐拉法變形的歐拉法中中 點(diǎn)點(diǎn) 方方 法法122222112121ccc ac b 36三階龍格三階龍格-庫塔方法庫塔方法三階龍格三階龍格-庫塔方法是用三個(gè)值庫塔方法是用三個(gè)值 k1, k2, k3 的加權(quán)平均來的加權(quán)平均來近似取代近似取代 k*111232213313231213122()(,)(,)(,)nnnnnnnnyyhkkkkf xykf xh yhkkf xh yhcccababkb hk 要使三階龍格要使三階龍格-庫塔方法具有三階精度，必須使其局部庫塔方法具有三階精度，必

32、須使其局部截?cái)嗾`差為截?cái)嗾`差為 O(h4)將將 k1, k2, k3 代入代入 yn 1 的表達(dá)式中，在的表達(dá)式中，在 (xn, , yn) 處用二元處用二元泰勒公式展開，與泰勒公式展開，與 y(xn 1) 在在 xn 處的泰勒展開式比較處的泰勒展開式比較37三階龍格三階龍格-庫塔方法（續(xù)）庫塔方法（續(xù)）類似二階龍格類似二階龍格-庫塔方法的推導(dǎo)過程，庫塔方法的推導(dǎo)過程，8 個(gè)待定系數(shù)個(gè)待定系數(shù) c1, c2, c3, a2, a3, b21, b31, b32 應(yīng)滿足：應(yīng)滿足： 123221331322233222233232211 21 31 6cccababbc ac ac ac ac b

33、 a 8 個(gè)未知參數(shù)，個(gè)未知參數(shù)，6 個(gè)個(gè)方程，有無窮多組解方程，有無窮多組解11231213121416661122()(,)(,)()2,nnnnnnnnyyhkkkkf xykf xh yhkkf xh yhkhk 庫塔公式庫塔公式38四階龍格四階龍格-庫塔方法庫塔方法類似可以推出四階龍格類似可以推出四階龍格-庫塔公式，常用的有：庫塔公式，常用的有：1123122166661122141213122243()(,)(,)(,)(,)nnnnnnnnnnyyhkkkkkf xykf xh yhkkf xh yhkkf xh yhk 標(biāo)準(zhǔn)四階龍格標(biāo)準(zhǔn)四階龍格- -庫塔公式庫塔公式39四階龍格

34、四階龍格-庫塔方法（續(xù)）庫塔方法（續(xù)）222211666611222 12212222222112341213124232()(,)(,)(,)(,)nnnnnnnnnnyyhkkkkkf xykf xh yhkkf xh yhkhkkf xh yhkhk 吉爾（吉爾（Gill）公式）公式4 階以上階以上龍格龍格-庫塔公式的計(jì)算量太大，并且精度不庫塔公式的計(jì)算量太大，并且精度不一定提高，有時(shí)反而會(huì)降低，因此實(shí)際應(yīng)用中一般一定提高，有時(shí)反而會(huì)降低，因此實(shí)際應(yīng)用中一般選用四階龍格選用四階龍格-庫塔已足可滿足精度要求。庫塔已足可滿足精度要求。40用經(jīng)典四階龍格用經(jīng)典四階龍格- -庫塔庫塔方法求解前例

35、的初值問題，并方法求解前例的初值問題，并與改進(jìn)與改進(jìn) 歐拉歐拉 法、梯形法在法、梯形法在 x5 0.5 處比較其誤差大小處比較其誤差大小解解：采用經(jīng)典四階龍格：采用經(jīng)典四階龍格- -庫塔庫塔公式：公式：10000(,)10kf xyyx 10.1,0.5(0)1yyxhxy 計(jì)計(jì)算算到到112220011122010(,)()()10.05kf xh yhkyhkxh 112230021122020(,)()()10.0475kf xh yhkyhkxh 4003030(,)()()10.09525kf xh yhkyhkxh 411221666610123412216666()10.1(00

36、.050.04750.09525)1.0048375yyhkkkk 于于是是：23451.01873091.040818421.070320291.1065309yyyy 同同理理可可計(jì)計(jì)算算：四階四階R-K方法方法的精度比二階的精度比二階方法高得多方法高得多10.1,0.5(0)1yyxhxy 計(jì)計(jì)算算到到0.55( )e()0.5e1.106530660 xy xxy x 精確解為：精確解為：R-K方法的誤差：方法的誤差：755()2.4 10yy x 改進(jìn)歐拉法的誤差：改進(jìn)歐拉法的誤差：45.5 10 梯形法的誤差：梯形法的誤差：42.5 10 42變步長的龍格變步長的龍格-庫塔方法庫塔

37、方法( )5111()hnnnnTy xyc h 設(shè)設(shè) y (xn) 在在 xn 處的值處的值 yn y (xn)，當(dāng)，當(dāng) xn+1 xn+ h 時(shí)時(shí) y (xn 1) 的近似值為的近似值為 ，由于四階，由于四階 R-K 方法的精度方法的精度為為 4 階，故局部截?cái)嗾`差為：階，故局部截?cái)嗾`差為：( )1hny 用四階用四階R-K方法求解初值問題精度較高，但要從理論方法求解初值問題精度較高，但要從理論上給出誤差上給出誤差 | | y (xn) yn | | 的估計(jì)式則比較困難；那么的估計(jì)式則比較困難；那么應(yīng)如何判斷計(jì)算結(jié)果的精度以及如何選擇合適的步應(yīng)如何判斷計(jì)算結(jié)果的精度以及如何選擇合適的步長長

38、 h？通常是通過不同步長在計(jì)算機(jī)上的計(jì)算結(jié)果進(jìn)行近通常是通過不同步長在計(jì)算機(jī)上的計(jì)算結(jié)果進(jìn)行近似估計(jì)。似估計(jì)。43若以若以 h/2 為步長，從為步長，從 xn 出發(fā)，經(jīng)過兩步計(jì)算，得到出發(fā)，經(jīng)過兩步計(jì)算，得到y(tǒng)(xn+1) 的近似值的近似值(2)1hny 變步長的龍格變步長的龍格-庫塔方法（續(xù)）庫塔方法（續(xù)）以上每步的截?cái)嗾`差約為以上每步的截?cái)嗾`差約為 cn(h/2)5，于是兩步的局部，于是兩步的局部截?cái)嗾`差為：截?cái)嗾`差為：(2)5111()2(2)hnnnnTy xyc h(2)511( )511()2(2)1()16hnnnhnnny xychy xyc h 于是：于是：整理得：整理得：(

39、2)(2)( )11111()15hhhnnnny xyyy44變步長的龍格變步長的龍格-庫塔方法（續(xù)）庫塔方法（續(xù)）記：記： ，給定的精度要求為，給定的精度要求為 e e(2)( )11115hhnnyy (2)1hny e e，反復(fù)將步長折半計(jì)算，直至，反復(fù)將步長折半計(jì)算，直至 e e，取最終得，取最終得到的到的 作為作為 y(xn+1) 的近似值。的近似值。 e e，反復(fù)將步長加倍計(jì)算，直至，反復(fù)將步長加倍計(jì)算，直至 e e，再將步長，再將步長折半一次計(jì)算，最終得到符合精度要求的折半一次計(jì)算，最終得到符合精度要求的 y(xn+1) 的的近似值。近似值。45單步法的收斂性單步法的收斂性顯式

40、單步法可統(tǒng)一寫成：顯式單步法可統(tǒng)一寫成：1(, )nnnnyyhxy h 增量函數(shù)，僅依賴于函數(shù)增量函數(shù)，僅依賴于函數(shù) f，且僅僅是且僅僅是 xn, yn, h 的函數(shù)的函數(shù)求求 y y(x)求求 y(xn)，xn x0 nh離散化離散化求求 yn y(xn)某種數(shù)值方法某種數(shù)值方法h 0 時(shí)，近似解時(shí)，近似解是否收斂到精確解是否收斂到精確解，它應(yīng)當(dāng)，它應(yīng)當(dāng)是一個(gè)固定節(jié)點(diǎn)，因是一個(gè)固定節(jié)點(diǎn)，因此此 h 0 時(shí)應(yīng)同時(shí)附時(shí)應(yīng)同時(shí)附帶帶 n 0nxxnh ( , , )x y h 46單步法的收斂性（續(xù)）單步法的收斂性（續(xù)）對(duì)于對(duì)于 p 階的常微分方程數(shù)值算法，當(dāng)階的常微分方程數(shù)值算法，當(dāng) h 0，

41、 n 時(shí)，是否時(shí)，是否 yn+1 y(xn+1)？p 階算法的局部截?cái)嗾`差為：階算法的局部截?cái)嗾`差為：111()()pnny xyO h 顯然：顯然：110,lim()0nnhny xy局部截?cái)嗾`差的前提假設(shè)是：局部截?cái)嗾`差的前提假設(shè)是： ()nnyy x 局部截?cái)嗾`差局部截?cái)嗾`差 0 并不能保證算法收斂并不能保證算法收斂47單步法的收斂性（續(xù)）單步法的收斂性（續(xù)）定義定義：若求解某初值問題的單步數(shù)值法，對(duì)于固定的：若求解某初值問題的單步數(shù)值法，對(duì)于固定的 當(dāng)當(dāng) h 0 且且 n 時(shí)，它的近似時(shí)，它的近似 解趨向于精確解解趨向于精確解 y(xn)，即：，即：0nxxnh 0,()limnnhn

42、yy x 則稱該單步法是收斂的。則稱該單步法是收斂的。定義定義：稱：稱 y(xn) yn 為單步法的近似解為單步法的近似解 yn 的的整體截?cái)嗾w截?cái)?誤差誤差。單步法收斂單步法收斂h 0 且且 n 時(shí)，時(shí)，yn 的整體截?cái)嗾`差的整體截?cái)嗾`差 0 48單步法的收斂性（續(xù)）單步法的收斂性（續(xù)）收斂性定理收斂性定理若某單步法滿足以上條件，則該方法是收斂的若某單步法滿足以上條件，則該方法是收斂的則該單步法的整體截?cái)嗾`差為：則該單步法的整體截?cái)嗾`差為：若單步法若單步法 具有具有 p 階精度，階精度，且增量函數(shù)且增量函數(shù) 關(guān)于關(guān)于 y 滿足：滿足：1(, )nnnnyyhxy h ( , , )x y

43、h ()()pnny xyO h ( , , )( , , )x y hx y hL yy Lipschitz 條件：條件： 初值初值 y0 是準(zhǔn)確的是準(zhǔn)確的49假設(shè)在前一步假設(shè)在前一步 yn 準(zhǔn)確的前提下求得的近似值為：準(zhǔn)確的前提下求得的近似值為：( , , )( , , )x y hx y hL yy 1(, )nnnnyyhxy h 1,()()nnnnyhxhy xy x 算法精度為算法精度為 p 階，局部截?cái)嗾`差：階，局部截?cái)嗾`差：111()pnny xych 111111() ()nnnnnny xyy xyyy 1111()nnnny xyyy 1 (), (), )(, )pn

44、nnnnnchy xyhxy xhhxy h 1(), (), )(, )pnnnnnnchy xyhxy xhxyh 1()()pnnnnchy xyhL y xy 1(1)()pnnhLy xych ne e1ne e 5011(1)pnnhLch e ee e 11(1)pnnhLch e ee e 11(1) (1)ppnhLhLchch e e 2 221(1)(1)1pnhLhLch e e 2 22113(1)(1)(1)1ppnhLhLchhLch e e 3213(1)(1)(1)1pnhLhLhLch e e 110(1)(1)(1)1nnphLhLhLch e e 101

45、(1)1(1)(1)1nnphLhLchhL e e 0(1)(1)1pnnchhLhLL e e 510(1)(1)1pnnnchhLhLL e ee e 2ee112xxxx (1)ennxx 即：即： (1)eeenhLnhLTLnhL 0(1)e1pTLnnchhLL e ee e 若初值是準(zhǔn)確的，則若初值是準(zhǔn)確的，則 e e 0 0 ，從而整體截?cái)嗾`差為：，從而整體截?cái)嗾`差為： e1pTLnchL e e 0nxxnhTy e x 為單調(diào)增函數(shù)，當(dāng)為單調(diào)增函數(shù)，當(dāng) 時(shí)時(shí)當(dāng)當(dāng) h 0 且且 n 時(shí)時(shí)，則，則0ne e()pO h52單步法的穩(wěn)定性單步法的穩(wěn)定性在討論單步法收斂性時(shí)一般認(rèn)

46、為數(shù)值方法本身的計(jì)在討論單步法收斂性時(shí)一般認(rèn)為數(shù)值方法本身的計(jì)算過程是準(zhǔn)確的，實(shí)際上并非如此：算過程是準(zhǔn)確的，實(shí)際上并非如此：初始值初始值 y0 有誤差有誤差 d d y0 y(x0)后續(xù)的每一步計(jì)算均有舍入誤差后續(xù)的每一步計(jì)算均有舍入誤差這些初始和舍入誤差在計(jì)算過程的傳播中是逐步衰這些初始和舍入誤差在計(jì)算過程的傳播中是逐步衰減的還是惡性增長就是數(shù)值方法的穩(wěn)定性問題減的還是惡性增長就是數(shù)值方法的穩(wěn)定性問題53定義定義：若一種數(shù)值方法在節(jié)點(diǎn)：若一種數(shù)值方法在節(jié)點(diǎn) xn 處的數(shù)值解處的數(shù)值解 yn 的擾的擾動(dòng)動(dòng) ，而在以后各節(jié)點(diǎn)，而在以后各節(jié)點(diǎn) ym (m n) 上產(chǎn)生的上產(chǎn)生的擾動(dòng)為擾動(dòng)為 ，如

47、果：，如果：md d0nd d 單步法的穩(wěn)定性（續(xù)）單步法的穩(wěn)定性（續(xù)）(1,2,)mnmnnd dd d 定義定義：設(shè)在節(jié)點(diǎn)：設(shè)在節(jié)點(diǎn) xn 處用數(shù)值算法得到的理想數(shù)值解為處用數(shù)值算法得到的理想數(shù)值解為 yn，而實(shí)際計(jì)算得到的近似解為，而實(shí)際計(jì)算得到的近似解為 ，稱差值：，稱差值：ny nnnyyd d 為第為第 n 步的數(shù)值解的步的數(shù)值解的擾動(dòng)擾動(dòng)。則稱該數(shù)值方法是則稱該數(shù)值方法是穩(wěn)定穩(wěn)定的。的。54單步法的穩(wěn)定性（續(xù)）單步法的穩(wěn)定性（續(xù)）歐拉法：歐拉法：1(,)nnnnyyhf xy 由于函數(shù)由于函數(shù) f (x, y) 的多樣性，數(shù)值穩(wěn)定性的分析相當(dāng)復(fù)的多樣性，數(shù)值穩(wěn)定性的分析相當(dāng)復(fù)雜，

48、通常只研究模型方程雜，通常只研究模型方程(0)yy 考察模型方程：考察模型方程：(0)yy 1(1)nnyhy 即：即：假設(shè)在節(jié)點(diǎn)值假設(shè)在節(jié)點(diǎn)值 yn 上有擾動(dòng)上有擾動(dòng) d dn，在節(jié)點(diǎn)值，在節(jié)點(diǎn)值 yn 1 上有上有擾動(dòng)擾動(dòng) d dn 1，且，且 d dn 1 僅由僅由 d dn 引起（即：計(jì)算過程中引起（即：計(jì)算過程中不再引起新的誤差）不再引起新的誤差）55111(1)(1)(1)()(1)nnnnnnnnyyhyhyhyyhd d d d 歐拉法穩(wěn)定歐拉法穩(wěn)定1nnd dd d 11h 即：即：111h 歐拉法穩(wěn)定的條件：歐拉法穩(wěn)定的條件：20h 0 1(1)nnyhy 針對(duì)模型方程：針

49、對(duì)模型方程：的顯式歐拉法：的顯式歐拉法：1(,)nnnnyyhf xy 化簡(jiǎn)得：化簡(jiǎn)得：20h 56隱式歐拉法：隱式歐拉法：111(,)nnnnyyhf xy 111111nnnnnnyyyyhhhd dd d 考察模型方程：考察模型方程：(0)yy 即：即：11nnnyyh y 化簡(jiǎn)為：化簡(jiǎn)為：11nnyyh 假設(shè)假設(shè) yn 上有擾動(dòng)上有擾動(dòng) ，則，則 yn 1 的擾動(dòng)為：的擾動(dòng)為：nnnyyd d 隱式歐拉法穩(wěn)定隱式歐拉法穩(wěn)定1nnd dd d 111h 0 0h ，上式均成立，所以：，上式均成立，所以：隱式歐拉法穩(wěn)定是恒穩(wěn)定的隱式歐拉法穩(wěn)定是恒穩(wěn)定的57一階常微分方程組一階常微分方程組1

50、1122212121012020( ,)( ,).( ,)(),(),()mmmmmmmyfx yyyyfx yyyyfx yyyy xs yxsyxs 111222( ,)( )( )( ,)( )( ,)( )( ,)TTTTmmmfx Yy xsyxsfx YY xF x YSyxsfx Y 0( )( ,)()TYxF x YY xS 5811121012212202( ,)()( ,)()yfx yyy xsyfx yyyxs 1(,)TnnnnYYhF x Y 111111222212(,)(,)nnnnnnnnnnyyfxyyhyyfxyy 111222( )( ,)( )( ,)( )( ,)TTTy xfx YsY xF x YSyxsfx Y 11(,)TnnnnYYF x Y111111111112222112(,)(,)nnnnnnnnnnyyfxyyhyyfxyy 111111111121112222122112