常微分方程數(shù)值解59382PPT學(xué)習(xí)教案

1、會計(jì)學(xué)1 常微分方程數(shù)值解常微分方程數(shù)值解59382 /* Eulers Method */ 歐拉公式：歐拉公式： x0 x1 向前差商近似導(dǎo)數(shù)向前差商近似導(dǎo)數(shù) h xyxy xy )()( )( 01 0 ),()()()( 000001 yxfhyxyhxyxy 1 y 記為記為 )1,., 0(),( 1 niyxfhyy iiii 定義定義在假設(shè)在假設(shè) yi = y(xi)，即第，即第 i 步計(jì)算是精確的前提下，考步計(jì)算是精確的前提下，考 慮的截?cái)嗾`差慮的截?cái)嗾`差 Ri = y(xi+1) yi+1 稱為稱為局部截?cái)嗾`差局部截?cái)嗾`差 /* local truncation error

2、*/。 定義定義若某算法的局部截?cái)嗾`差為若某算法的局部截?cái)嗾`差為O(hp+1)，則稱該算法有，則稱該算法有p 階精度。階精度。 歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差： ),()()()()()( 3 211 2 iiii h iiiii yxhfyhOxyxyhxyyxyR )()( 3 2 2 hOxy i h 歐拉法具有歐拉法具有 1 階精度。階精度。 Ri 的的主項(xiàng)主項(xiàng) /* leading term */ 亦稱為亦稱為歐拉折線法歐拉折線法 /* Eulers polygonal arc method*/ 第1頁/共19頁 解：解：h=0.2 , xi=1+ih， 精確解為：精確解

3、為：y=x2-2x y1= y0+hf(x0, y0)=-1 y2= y1+hf(x1, y1)=- 0.9333 y3= y2+hf(x2, y2)=-0.8 y4= y3+hf(x3, y3)=-0.6 y5= y4+hf(x4, y4)=-0.3333 y6= y5+hf(x5, y5)=0 dy2y 2 dxx y 11 可以看出誤差隨著計(jì)算在積累。可以看出誤差隨著計(jì)算在積累。 xky(xk)ykek 1.2- 0.96-10.04 1.4-0.84 - 0.9333 0.0933 1.6-0.64-0.80.16 1.8-0.36-0.60.24 2.00 - 0.3333 0.33

4、33 2.20.4400.44 步長步長h=0.2 第2頁/共19頁 1 Eulers Method 隱式歐拉法隱式歐拉法 /* implicit Euler method */ 向后差商近似導(dǎo)數(shù)向后差商近似導(dǎo)數(shù) h xyxy xy )()( )( 01 1 x0 x1 )(,()( 1101 xyxfhyxy )1,., 0(),( 111 niyxfhyy iiii 由于未知數(shù)由于未知數(shù) yi+1 同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故 稱為稱為隱式隱式 /* implicit */ 歐拉公式，而前者稱為歐拉公式，而前者稱為顯式顯式 /* explic

5、it */ 歐拉公式。歐拉公式。 一般先用顯式計(jì)算一個(gè)初值，再一般先用顯式計(jì)算一個(gè)初值，再迭代迭代求解。求解。 隱式隱式歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差： 11) ( iii yxyR)()( 3 2 2 hOxy i h 即隱式歐拉公式具有即隱式歐拉公式具有 1 階精度。階精度。 第3頁/共19頁 1 Eulers Method /* trapezoid formula */ 顯顯/隱式兩種算法的隱式兩種算法的平均平均 )1,., 0(),(),( 2 111 niyxfyxf h yy iiiiii 注：注：的確有局部截?cái)嗾`差的確有局部截?cái)嗾`差 ， 即梯形公式具有即梯形公式具有2

6、 階精度，比歐拉方法有了進(jìn)階精度，比歐拉方法有了進(jìn) 步。但注意到該公式是步。但注意到該公式是隱式隱式公式，計(jì)算時(shí)不得公式，計(jì)算時(shí)不得 不用到迭代法。不用到迭代法。 )()( 3 11 hOyxyR iii 中點(diǎn)歐拉公式中點(diǎn)歐拉公式 /* midpoint formula */ 中心差商近似導(dǎo)數(shù)中心差商近似導(dǎo)數(shù) h xyxy xy 2 )()( )( 02 1 x0 x2x1 )(,(2)()( 1102 xyxfhxyxy 1,., 1),(2 11 niyxfhyy iiii 假設(shè)假設(shè) ，則可以導(dǎo)出，則可以導(dǎo)出 即中點(diǎn)公式具有即中點(diǎn)公式具有 2 階精度。階精度。 )(),( 11iiii x

7、yyxyy )()( 3 11 hOyxyR iii 需要需要2個(gè)初值個(gè)初值 y0和和 y1來啟動遞推來啟動遞推 過程，這樣的算法稱為過程，這樣的算法稱為雙步法雙步法 /* double-step method */，而前面的三種算法都是，而前面的三種算法都是單步法單步法 /* single-step method */。 第4頁/共19頁 方方 法法 1 Eulers Method 顯式歐拉顯式歐拉 隱式歐拉隱式歐拉 梯形公式梯形公式 中點(diǎn)公式中點(diǎn)公式 簡單簡單精度低精度低 穩(wěn)定性最好穩(wěn)定性最好精度低精度低, 計(jì)算量大計(jì)算量大 精度提高精度提高計(jì)算量大計(jì)算量大 精度提高精度提高, 顯式顯式

8、多一個(gè)初值多一個(gè)初值, 可能影響精度可能影響精度 Cant you give me a formula with all the advantages yet without any of the disadvantages? Do you think it possible? Well, call me greedy OK, lets make it possible. 第5頁/共19頁 /* modified Eulers method */ Step 1: 先用先用顯式顯式歐拉公式作歐拉公式作預(yù)測預(yù)測，算出，算出 ),( 1iiii yxfhyy Step 2: 再將再將 代入代入隱式隱

9、式梯形公式的右邊作梯形公式的右邊作校正校正，得到，得到 1 i y ),(),( 2 111 iiiiii yxfyxf h yy 注：注：此法亦稱為此法亦稱為預(yù)測預(yù)測-校正法校正法 /* predictor-corrector method */。 可以證明該算法具有可以證明該算法具有 2 階精度，同時(shí)可以看到它是個(gè)階精度，同時(shí)可以看到它是個(gè) 單步單步遞推格式，比隱式公式的迭代求解過程遞推格式，比隱式公式的迭代求解過程簡單簡單。后。后 面將看到，它的面將看到，它的穩(wěn)定性高穩(wěn)定性高于顯式歐拉法。于顯式歐拉法。 ) 1,., 0(),(,),( 2 11 niyxfhyxfyxf h yy ii

10、iiiiii 1 Eulers Method 第6頁/共19頁 解：梯形公式為：解：梯形公式為： , ( ) dy2x yy 01 dxy kk 1 k 1kkk 1 kk 1 h2x2x yyyy 2yy ()() 2k k 1 k 1kk 1 k hhhx y1y1yhx0 22y 可以預(yù)測校正方法來求：可以預(yù)測校正方法來求： () () () k pkk k k 1 ckp p k 1pc 2x yyh y y 2x yyh y y 1 yyy 2 解解yk+1出來比較困難，遇到的是一個(gè)二次方程，出來比較困難，遇到的是一個(gè)二次方程， 即：即： ),(,),( 2 11kkkkkkkk y

11、xfhyxfyxf h yy 第7頁/共19頁 xkyky(xk)ek 0.11.09596911.09544510.000464 0.21.18409661.18321600.000881 0.31.26620141.26491110.001290 0.41.34336021.34164080.001719 0.51.41640191.41421360.002188 0.61.47595561.48323970.002716 0.71.55251411.54919330.003321 0.81.61647481.61245150.004023 0.91.67816641.67332010.0

12、04746 1.01.73786741.73205080.005817 ),(,),( 2 11kkkkkkkk yxfhyxfyxf h yy 第8頁/共19頁 /* Runge-Kutta Method */ 建立高精度的單步遞推格式。建立高精度的單步遞推格式。 單步遞推法的單步遞推法的基本思想基本思想是從是從 ( 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)的歐拉法，可以將其改寫為：

13、),( ),( 2 1 2 1 12 1 211 hKyhxfK yxfK KKhyy ii ii ii 斜率斜率 一定取一定取K1 K2 的的平均值平均值嗎？嗎？ 步長一定是一個(gè)步長一定是一個(gè)h 嗎嗎 ？ 第9頁/共19頁 2 Runge-Kutta Method 首先希望能確定系數(shù)首先希望能確定系數(shù) 1、 2、p，使得到的算法格式有，使得到的算法格式有2階階 精度，即在精度，即在 的前提假設(shè)下，使得的前提假設(shè)下，使得 )( ii xyy )()( 3 11 hOyxyR iii Step 1: 將將 K2 在在 ( xi , yi ) 點(diǎn)作點(diǎn)作 Taylor 展開展開 )(),(),(),

14、( ),( 2 1 12 hOyxfphKyxphfyxf phKyphxfK iiyiixii ii )()()( 2 hOxyphxy ii 將改進(jìn)歐拉法推廣為：將改進(jìn)歐拉法推廣為： ),( ),( 12 1 22111 phKyphxfK yxfK KKhyy ii ii ii ),(),(),( ),(),( ),()( yxfyxfyxf dx dy yxfyxf yxf dx d xy yx yx Step 2: 將將 K2 代入第代入第1式，得到式，得到 )()()()( )()()()( 32 221 2 211 hOxyphxyhy hOxyphxyxyhyy iii iii

15、ii 第10頁/共19頁 2 Runge-Kutta Method Step 3: 將將 yi+1 與與 y( xi+1 ) 在在 xi 點(diǎn)的點(diǎn)的泰勒泰勒展開作比較展開作比較 )()()()( 32 2211 hOxyphxyhyy iiii )()( 2 )()()( 3 2 1 hOxy h xyhxyxy iiii 要求要求 ，則必須有，則必須有 ： )()( 3 11 hOyxyR iii 2 1 ,1 221 p 這里有這里有 個(gè)未知個(gè)未知 數(shù)，數(shù)， 個(gè)方程。個(gè)方程。 3 2 存在存在無窮多個(gè)解無窮多個(gè)解。所有滿足上式的格式統(tǒng)稱為。所有滿足上式的格式統(tǒng)稱為2階龍格階龍格 - 庫塔格式

16、庫塔格式。 2 1 , 1 21 p注意到，注意到， 就是改進(jìn)的歐拉法。就是改進(jìn)的歐拉法。 Q: 為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？ 第11頁/共19頁 其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均為均為 待定系數(shù)，確定這些系待定系數(shù)，確定這些系 數(shù)的步驟與前面相似。數(shù)的步驟與前面相似。 2 Runge-Kutta Method ).,( . ),( ),( ),( . 112211 23213133 12122 1 22111 mm mmmmi

17、m ii ii ii mmii hKhKhKyihxfK hKhKyhxfK hKyhxfK yxfK KKKhyy 最常用為四級最常用為四級4階階經(jīng)典龍格經(jīng)典龍格-庫塔法庫塔法 /* Classical Runge-Kutta Method */ ： ),( ),( ),( ),( )22( 34 2223 1222 1 432161 hKyhxfK KyxfK KyxfK yxfK KKKKyy ii h i h i h i h i ii h ii 第12頁/共19頁 0021737. 2 2 2 2 2 2 , 2 20 0 0 20 2003 K h y h x h x K h y K

18、 h y h xfK 007603. 2, 30 0 0 30 3004 hKy hx hx hKy hKyhxfK 2002714. 1 32 6 432101 KKKK h yy 11y y x x y dx dy 1ln2xxxy 其精確解為：其精確解為： 用經(jīng)典龍格庫塔方法用經(jīng)典龍格庫塔方法，h=0.1 1 . 0, 1, 1 00 hxy(1) 求求y1，此時(shí)，此時(shí) 2, 0 0 0 0 001 y x x y yxfK 0021645. 2 2 2 2 2 2 , 2 10 0 0 10 1002 K h y h x h x K h y K h y h xfK 解：解： 20027

19、11. 11ln2 111 xxxy 第13頁/共19頁 k x k y k xy k e 7 103 7 105 7 106 1.11.200271 4 1.2002711 1.21.401815 8 1.4018153 1.31.605239 3 1.6052387 1.41.810793 6 1.7107930 1.52.018562 8 2.0185621 1.62.228547 0 2.2285463 1.72.440703 8 2.4407030 1.82.654969 3 2.6549685 1.92.871269 7 2.8712689 2.03.089527 9 3.0895

20、271 7 106 7 107 7 107 7 108 7 108 7 108 7 108 第14頁/共19頁 1 22211 33311322 1112233 (,) (,) (,() () kk kk kk kk Kf xy Kf xa h yhK Kf xa h yhKK yyhKKK 共有八個(gè)參數(shù)：共有八個(gè)參數(shù)： 1, 2, 3, 1, 2, 21, 31, 32， 1、三階龍格庫塔公式格式如下：、三階龍格庫塔公式格式如下： v 將將 f 在在 xk 點(diǎn)點(diǎn)泰勒展開泰勒展開： 1 , kkk Kf xyyx 2222223 2221222121 11 () 22 xyxxxyyy Kfa hfhffa h fah ffh f fO h 222 3331132233311322 2 3 311322 1 () 2 1 2 xyxxxy yy KfhfKKhffh fKKh f KKfh 第15頁/共19頁 將將K1,K2按展開式代入按展開式