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

1、三常微分方程數(shù)值解法2常微分方程數(shù)值解法引言（常微分方程數(shù)值解法概述）引言（常微分方程數(shù)值解法概述）顯式歐拉法、隱式歐拉法、二步歐拉法顯式歐拉法、隱式歐拉法、二步歐拉法局部截斷誤差與精度局部截斷誤差與精度改進的歐拉方法改進的歐拉方法龍格龍格-庫塔方法庫塔方法收斂性與穩(wěn)定性簡述收斂性與穩(wě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，使得對于，使得對于 R 域內(nèi)的任意兩值域內(nèi)的任意兩值 y1, y2，下，下列不等式成立：列不等式成立： 1221|( ,)( ,)|f x yf x yL yy則上述初值問題的連續(xù)可微的解則上述初值問題的連續(xù)可微的解 y(x) 存在并且唯一。存在并且唯一。4引言（續(xù)）引言（續(xù)）實際生產(chǎn)與科研中，除少數(shù)簡單情況能獲得初值問題實際生產(chǎn)與科研中，除少數(shù)簡單情況能獲得初值問題的初等解（用初等函數(shù)表示的解）外，絕大多數(shù)情況的初等解（用初等函數(shù)表示的解）外，絕大多數(shù)情況下是求不出初等解的。下

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

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

5、函數(shù)值數(shù)據(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ū)間寬度與右端點處的函數(shù)值乘積，即：度與右端點處的函數(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 需要前兩步需要前兩步的計算結(jié)果的計算結(jié)果中矩形公式中矩形公式11梯形公式歐拉法梯形公式歐拉法在數(shù)值積分法中，如果用梯形公式近似計算在數(shù)值積分法中，如果用梯形公式近似計算 f (x, y) 在區(qū)間在區(qū)間 xn, xn+1 上的積分，即：上的積分，即：用近似值代替精確值可得梯形公式歐拉法：用近似值代替精確值可得梯形公式歐拉法：111 (,)(,)2nnnnnnhyyf xyf xy上式右端出現(xiàn)了未知項，可見梯形法是隱式歐拉法上式右端出現(xiàn)了未知項，可見梯形法是隱式歐拉法的一種；實際上，梯形公式歐拉法是顯式歐拉法與的一種；實際上，梯形公式歐

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，計算到，計算到 x 0.5，并與精確解進行比較，并與精確解進行比較解：由已知條件可得：解：由已知條件可得：h 0.1，x0 0， y0 1， f (x, y) y x 1顯式歐拉法：顯式歐拉法：1(1)(1)0.90.10.1nnnn

10、nnnnyyhyxh yhxhyx 13例：（續(xù)）例：（續(xù)）隱式歐拉法：隱式歐拉法：111(1)nnnnyyhyx化簡得：化簡得：11(1)(1)(1)nnnnnh yyh xyh xh11(1)(0.10.11) 1.1(1)nnnnnyyh xhyxh 111(1)(1)2nnnnnnhyyyxyx 梯形公式歐拉法：梯形公式歐拉法：111(2)1(2)(2)(2)nnnnhyhyhxx11(2)(2) (2)(1.90.20.21) 2.1nnnnnnyh yh xxhyx14計算結(jié)果：計算結(jié)果：xn顯式法顯式法 yn隱式法隱式法 yn梯形法梯形法 yn精確解精確解 y (xn)0.011

11、110.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局部截斷誤差局部截斷誤差為了簡

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

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

14、開式的前幾項將將 y(xn 1) 在在 xn 點處用泰勒公式展開：點處用泰勒公式展開：211( )()()()()2!(,)nnnnnnyy xy xhy xhy xhxx 22111( )()()2!nnnyTy xyhO h 顯式歐拉法的局部截斷誤差：顯式歐拉法的局部截斷誤差：1(), ()()()nnnnnnyy xhf xy xy xhy x 1(,)nnnnyyhf xy 歐拉法歐拉法()nnyy x 令令：1 階精度階精度18補充：二元函數(shù)微分中值定理補充：二元函數(shù)微分中值定理1100(,)(,)f xyf xy 1001010010()()()()xyxxfxxxyyfyyy 0

15、1 0000(,)(,)f xh ykf xy00()()xyhfxhkfyk000000(,)(,)()()xyf xh ykf xyhfxhkfyk0000(,)(,)f xyhkf xh ykxy191111111()(,),(, )()nnnynnnnyf xyf xxxfxyy 111()()( )(, )nnnynnyy xh y xO hfxT 21()()()()nnny xy xhy xO h y(xn 1) 在在 xn 點處展開：點處展開：211111()(, )()nnnynnTy xyhfxTO h 111(,)nnnnyyhf xy隱式歐拉法：隱式歐拉法：1 階精度階

16、精度()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 點處用泰勒公式展開：點處用泰勒公式展開：231231()()()()()2!()()()()()()2!nnnnnnnnyxy xy xhy xhO hyxy xy xhy xhO h 111()nnnTy xy二

17、步歐拉法的局部截斷誤差：二步歐拉法的局部截斷誤差：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 點處展開：點處展開：231()()()()()

18、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各種歐拉法的比較各種歐拉法的比較方法方法精度精度評述評述顯式歐拉法顯式歐拉法1最簡單，精度低最簡單，精度低隱式歐拉法隱式歐拉法1不便計算，穩(wěn)定性好不便計算，穩(wěn)定性好二步歐拉法二步歐拉法2需要兩步初值，且第需要兩步初值，且第 2 個初值只能由其個初值只能由其它方法給出，

19、可能對后面的遞推精度有它方法給出，可能對后面的遞推精度有影響影響梯形公式梯形公式歐歐 拉拉 法法2精度有所提高，但為隱式，需要迭代求精度有所提高，但為隱式，需要迭代求解，計算量大解，計算量大23改進的歐拉法改進的歐拉法從上述例子可以看到，梯形法由于具有二階精度，其從上述例子可以看到，梯形法由于具有二階精度，其局部截斷誤差比顯式歐拉法和隱式歐拉法小，但梯形局部截斷誤差比顯式歐拉法和隱式歐拉法小，但梯形法實質(zhì)上是一種隱式算法法實質(zhì)上是一種隱式算法顯式歐拉法是一個顯式算法，雖然計算量較小，但是顯式歐拉法是一個顯式算法，雖然計算量較小，但是精度不高精度不高綜合兩種方法的長處，可以先用顯式歐拉法求出綜合

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

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

22、nhyyxyh f xyx (1)(1)12nnnnnnnhyyxyhyxxh( 2)(2)122nnhhhhyxh 0.9050.0950.1nnyx26計算結(jié)果計算結(jié)果xn改進歐改進歐拉法拉法 yn精確解精確解 y (xn) 誤誤 差差改進歐拉法改進歐拉法歐拉法歐拉法梯形法梯形法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.4 1.070802 1.0703204.

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

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

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

26、xn, xn 1 內(nèi)多預(yù)報幾個點的內(nèi)多預(yù)報幾個點的 ki 值，并用其加權(quán)平均值，并用其加權(quán)平均值作為值作為 k* 的近似值從而構(gòu)造出具有更高精度的計算公的近似值從而構(gòu)造出具有更高精度的計算公式，這就是式，這就是龍格龍格- -庫塔方法的基本思想。庫塔方法的基本思想。31二階龍格二階龍格-庫塔方法庫塔方法以以 k1 和和 k2 的加權(quán)平均來近似取代的加權(quán)平均來近似取代 k*12221121211()(,)(,)nnnnnncyyhkkkf xykf xh yhcakb 12221,ccab為為待待定定系系數(shù)數(shù)為分析局部截斷誤差，令為分析局部截斷誤差，令 yn y(xn)，由泰勒公式得：，由泰勒公式

27、得：231()()()()()()2!nnnnnhy xy xhy xhy xyxO h (), ()(,)nnnnny xf xy xf xy 123()()(,)(,)(,) (,)()2nnnnxnnynnnny xy xhf xyhfxyfxyf xyO h ()(,)(,)(,) (,)nnnxnnynnnnyxfxyfxyfxyf xyny 32補充：二元泰勒展開式補充：二元泰勒展開式00000020000(,)(,)(,)1(,)2!1(,)!nf xh ykf xyhkf xyxyhkf xyxyhkf xynxy 20000200(,)2(,)(,)xxxyyyh fxyhk

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

29、中可得：中可得：11122()nnyyh c kc k 2122223221() (,)(,)(,) (,)()nnnxnnynnnnyh ccf xyh c a fxyh c b fxyf xyO h122222213() (,)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()()(,)(,)(,) (,)()2nnnnxnnynnnny

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

31、,),22nnnnnnyyhkkf xyhhkfxyk 變形的歐拉法變形的歐拉法中中 點點 方方 法法122222112121ccc ac b36三階龍格三階龍格-庫塔方法庫塔方法三階龍格三階龍格-庫塔方法是用三個值庫塔方法是用三個值 k1, k2, k3 的加權(quán)平均來的加權(quán)平均來近似取代近似取代 k*111232213313231213122()(,)(,)(,)nnnnnnnnyyhkkkkf xykf xh yhkkf xh yhcccababkb hk 要使三階龍格要使三階龍格-庫塔方法具有三階精度，必須使其局部庫塔方法具有三階精度，必須使其局部截斷誤差為截斷誤差為 O(h4)將將 k

32、1, k2, k3 代入代入 yn 1 的表達式中，在的表達式中，在 (xn, , yn) 處用二元處用二元泰勒公式展開，與泰勒公式展開，與 y(xn 1) 在在 xn 處的泰勒展開式比較處的泰勒展開式比較37三階龍格三階龍格-庫塔方法（續(xù)）庫塔方法（續(xù)）類似二階龍格類似二階龍格-庫塔方法的推導(dǎo)過程，庫塔方法的推導(dǎo)過程，8 個待定系數(shù)個待定系數(shù) c1, c2, c3, a2, a3, b21, b31, b32 應(yīng)滿足：應(yīng)滿足： 123221331322233222233232211 21 31 6cccababbc ac ac ac ac b a 8 個未知參數(shù)，個未知參數(shù)，6 個個方程，有

33、無窮多組解方程，有無窮多組解11231213121416661122()(,)(,)()2,nnnnnnnnyyhkkkkf xykf xh yhkkf xh yhkhk 庫塔公式庫塔公式38四階龍格四階龍格-庫塔方法庫塔方法類似可以推出四階龍格類似可以推出四階龍格-庫塔公式，常用的有：庫塔公式，常用的有：1123122166661122141213122243()(,)(,)(,)(,)nnnnnnnnnnyyhkkkkkf xykf xh yhkkf xh yhkkf xh yhk 標準四階龍格標準四階龍格- -庫塔公式庫塔公式39四階龍格四階龍格-庫塔方法（續(xù)）庫塔方法（續(xù)）222211

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

35、進 歐拉歐拉 法、梯形法在法、梯形法在 x5 0.5 處比較其誤差大小處比較其誤差大小解解：采用經(jīng)典四階龍格：采用經(jīng)典四階龍格- -庫塔庫塔公式：公式：10000(,)10kf xyyx 10.1,0.5(0)1yyxhxy 計計算算到到112220011122010(,)()()10.05kf xh yhkyhkxh 112230021122020(,)()()10.0475kf xh yhkyhkxh 4003030(,)()()10.09525kf xh yhkyhkxh 411221666610123412216666()10.1(00.050.04750.09525)1.004837

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

37、設(shè) y (xn) 在在 xn 處的值處的值 yn y (xn)，當，當 xn+1 xn+ h 時時 y (xn 1) 的近似值為的近似值為 ，由于四階，由于四階 R-K 方法的精度方法的精度為為 4 階，故局部截斷誤差為：階，故局部截斷誤差為：( )1hny 用四階用四階R-K方法求解初值問題精度較高，但要從理論方法求解初值問題精度較高，但要從理論上給出誤差上給出誤差 | | y (xn) yn | | 的估計式則比較困難；那么的估計式則比較困難；那么應(yīng)如何判斷計算結(jié)果的精度以及如何選擇合適的步應(yīng)如何判斷計算結(jié)果的精度以及如何選擇合適的步長長 h？通常是通過不同步長在計算機上的計算結(jié)果進行近通

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

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

40、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 時，近似解時，近似解是否收斂到精確解是否收斂到精確解，它應(yīng)當，它應(yīng)當是一個固定節(jié)點，因是一個固定節(jié)點，因此此 h 0 時應(yīng)同時附時應(yīng)同時附帶帶 n 0nxxnh( , , )x y h 46單步法的收斂性（續(xù)）單步法的收斂性（續(xù)）對于對于 p 階的常微分方程數(shù)值算法，當階的常微分方程數(shù)值算法，當 h 0， n 時，是否時，是否 yn+1 y(xn+1)？p

41、 階算法的局部截斷誤差為：階算法的局部截斷誤差為：111()()pnny xyO h 顯然：顯然：110,lim()0nnhny xy局部截斷誤差的前提假設(shè)是：局部截斷誤差的前提假設(shè)是： ()nnyy x 局部截斷誤差局部截斷誤差 0 并不能保證算法收斂并不能保證算法收斂47單步法的收斂性（續(xù)）單步法的收斂性（續(xù)）定義定義：若求解某初值問題的單步數(shù)值法，對于固定的：若求解某初值問題的單步數(shù)值法，對于固定的 當當 h 0 且且 n 時，它的近似時，它的近似 解趨向于精確解解趨向于精確解 y(xn)，即：，即：0nxxnh0,()limnnhnyy x 則稱該單步法是收斂的。則稱該單步法是收斂的。

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

43、, )x y hx y hL yy Lipschitz 條件：條件： 初值初值 y0 是準確的是準確的49假設(shè)在前一步假設(shè)在前一步 yn 準確的前提下求得的近似值為：準確的前提下求得的近似值為：( , , )( , , )x y hx y hL yy 1(, )nnnnyyhxy h 1,()()nnnnyhxhy xy x 算法精度為算法精度為 p 階，局部截斷誤差：階，局部截斷誤差：111()pnny xych 111111() ()nnnnnny xyy xyyy1111()nnnny xyyy1 (), (), )(, )pnnnnnnchy xyhxy xhhxy h 1(), ()

44、, )(, )pnnnnnnchy xyhxy xhxy h 1()()pnnnnchy xyhL y xy 1(1)()pnnhLy xych ne e1ne e 5011(1)pnnhLch eeee 11(1)pnnhLch eeee 11(1) (1)ppnhLhLchche e 2 221(1)(1)1pnhLhLche e 2 22113(1)(1)(1)1ppnhLhLchhLche e 3213(1)(1)(1)1pnhLhLhLche e 110(1)(1)(1)1nnphLhLhLche e 101(1)1(1)(1)1nnphLhLchhL e e 0(1)(1)1pnn

45、chhLhLL e e510(1)(1)1pnnnchhLhLL eeee2ee112xxxx (1)ennxx即：即： (1)eeenhLnhLTLnhL 0(1)e1pTLnnchhLL eeee若初值是準確的，則若初值是準確的，則 e e 0 0 ，從而整體截斷誤差為：，從而整體截斷誤差為： e1pTLnchL e e0nxxnhTy e x 為單調(diào)增函數(shù)，當為單調(diào)增函數(shù)，當 時時當當 h 0 且且 n 時時，則，則0ne e()pO h52單步法的穩(wěn)定性單步法的穩(wěn)定性在討論單步法收斂性時一般認為數(shù)值方法本身的計在討論單步法收斂性時一般認為數(shù)值方法本身的計算過程是準確的，實際上并非如此：

46、算過程是準確的，實際上并非如此：初始值初始值 y0 有誤差有誤差 d d y0 y(x0)后續(xù)的每一步計算均有舍入誤差后續(xù)的每一步計算均有舍入誤差這些初始和舍入誤差在計算過程的傳播中是逐步衰這些初始和舍入誤差在計算過程的傳播中是逐步衰減的還是惡性增長就是數(shù)值方法的穩(wěn)定性問題減的還是惡性增長就是數(shù)值方法的穩(wěn)定性問題53定義定義：若一種數(shù)值方法在節(jié)點：若一種數(shù)值方法在節(jié)點 xn 處的數(shù)值解處的數(shù)值解 yn 的擾的擾動動 ，而在以后各節(jié)點，而在以后各節(jié)點 ym (m n) 上產(chǎn)生的上產(chǎn)生的擾動為擾動為 ，如果：，如果：md d0nd d 單步法的穩(wěn)定性（續(xù)）單步法的穩(wěn)定性（續(xù)）(1,2,)mnmnn

47、dddd定義定義：設(shè)在節(jié)點：設(shè)在節(jié)點 xn 處用數(shù)值算法得到的理想數(shù)值解為處用數(shù)值算法得到的理想數(shù)值解為 yn，而實際計算得到的近似解為，而實際計算得到的近似解為 ，稱差值：，稱差值：ny nnnyyd d 為第為第 n 步的數(shù)值解的步的數(shù)值解的擾動擾動。則稱該數(shù)值方法是則稱該數(shù)值方法是穩(wěn)定穩(wěn)定的。的。54單步法的穩(wěn)定性（續(xù)）單步法的穩(wěn)定性（續(xù)）歐拉法：歐拉法：1(,)nnnnyyhf xy 由于函數(shù)由于函數(shù) f (x, y) 的多樣性，數(shù)值穩(wěn)定性的分析相當復(fù)的多樣性，數(shù)值穩(wěn)定性的分析相當復(fù)雜，通常只研究模型方程雜，通常只研究模型方程(0)yy 考察模型方程：考察模型方程：(0)yy 1(1)

48、nnyhy 即：即：假設(shè)在節(jié)點值假設(shè)在節(jié)點值 yn 上有擾動上有擾動 d dn，在節(jié)點值，在節(jié)點值 yn 1 上有上有擾動擾動 d dn 1，且，且 d dn 1 僅由僅由 d dn 引起（即：計算過程中引起（即：計算過程中不再引起新的誤差）不再引起新的誤差）55111(1)(1)(1)()(1)nnnnnnnnyyhyhyhyyhd d d d 歐拉法穩(wěn)定歐拉法穩(wěn)定1nndddd 11h 即：即：111h 歐拉法穩(wěn)定的條件：歐拉法穩(wěn)定的條件：20h 0 1(1)nnyhy 針對模型方程：針對模型方程：的顯式歐拉法：的顯式歐拉法：1(,)nnnnyyhf xy 化簡得：化簡得：20h 56隱式

49、歐拉法：隱式歐拉法：111(,)nnnnyyhf xy111111nnnnnnyyyyhhhd dd d 考察模型方程：考察模型方程：(0)yy 即：即：11nnnyyh y 化簡為：化簡為：11nnyyh 假設(shè)假設(shè) yn 上有擾動上有擾動 ，則，則 yn 1 的擾動為：的擾動為：nnnyyd d 隱式歐拉法穩(wěn)定隱式歐拉法穩(wěn)定1nndddd 111h 0 0h ，上式均成立，所以：，上式均成立，所以：隱式歐拉法穩(wěn)定是恒穩(wěn)定的隱式歐拉法穩(wěn)定是恒穩(wěn)定的57一階常微分方程組一階常微分方程組11122212121012020( ,)( ,).( ,)(),(),()mmmmmmmyfx yyyyfx

50、yyyyfx yyyy xsyxsyxs 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 Y11(,)TnnnnYYF x Y111111111112222112(,)(,)nnnnnnnnnnyyfxyyhyyfxyy 111111111121112222122112(,