數(shù)值計(jì)算方法第六章

1、第六章第六章 常微分方程初值問題的常微分方程初值問題的 數(shù)值解法數(shù)值解法 6.1 歐拉方法歐拉方法 6.2 計(jì)算公式的誤差分析計(jì)算公式的誤差分析 6.3 龍格庫塔方法龍格庫塔方法 6.4 向一階方程組與高階方程的推廣向一階方程組與高階方程的推廣 問題的提出問題的提出 的的數(shù)數(shù)值值解解 求求初初值值問問題題)16(, )( ),( 0 00 bxx yxy yxfy 數(shù)值求解方法數(shù)值求解方法 的的近近似似值值 數(shù)數(shù)值值逐逐個(gè)個(gè)求求解解出出節(jié)節(jié)點(diǎn)點(diǎn)上上的的函函 生生成成節(jié)節(jié)點(diǎn)點(diǎn)作作等等距距分分割割對(duì)對(duì)區(qū)區(qū)間間 )(,),(),( )2( ), 1 , 0( ,)1( 21 0 0 0 N i xy

2、xyxy N xb hniihxx bx 6.1 歐拉方法歐拉方法 6.1.1 歐拉公式與改進(jìn)歐拉公式歐拉公式與改進(jìn)歐拉公式 算法：算法： N xb hniihxx bx i 0 0 0 ), 1 ,0( , 生生成成節(jié)節(jié)點(diǎn)點(diǎn)作作等等距距分分割割對(duì)對(duì)區(qū)區(qū)間間 處的近似值處的近似值在在 處的精確值處的精確值在在 ii ii xxyy xxyxy )( )()( 1 11 )(,()()( )(,()( ),(, 1 1 i i i i i i x x ii x x x x ii dxxyxfxyxy dxxyxfdxxy yxfyxx 即即 ：積積分分上上對(duì)對(duì)在在 )46()(,()()( 1

3、1 i i x x ii dxxyxfxyxy 化化為為積積分分方方程程問問題題從從而而把把微微分分方方程程問問題題轉(zhuǎn)轉(zhuǎn) 選擇不同的數(shù)值積分公式來求選擇不同的數(shù)值積分公式來求 近似值就得到初值問題的各種數(shù)值解法近似值就得到初值問題的各種數(shù)值解法 1 )(,( i i x x dxxyxf 1.歐拉公式歐拉公式 )(,()()( 1iiii xyxfhxyxy 得得，時(shí)時(shí)由由 當(dāng)當(dāng) )46( )(,()(,( 1 ii x x xyxfhdxxyxf i i ox y ),(yxfY 0 x 1 x 1 i xbxn i x h )26()1, 1 ,0( )( ),( 00 1 Ni xyy

4、yxfhyy iiii 由由此此可可建建立立計(jì)計(jì)算算公公式式 這稱為這稱為歐拉公式歐拉公式 2.后退歐拉公式后退歐拉公式 ox y ),(yxfY 0 x 1 x 1 i xbxn i x h 得得，時(shí)時(shí)由由 當(dāng)當(dāng) )46( )(,()(,( 11 1 ii x x xyxfhdxxyxf i i )36()1, 1 , 0( )( ),( 00 111 Ni xyy yxfhyy iiii 由由此此可可建建立立另另一一公公式式 )(,()()( 111 iiii xyxfhxyxy 這稱為這稱為后退歐拉公式后退歐拉公式 后退歐拉公式是一個(gè)隱式公式，通常采用迭代法求解。后退歐拉公式是一個(gè)隱式公

5、式，通常采用迭代法求解。 比比較較并并與與精精確確解解的的數(shù)數(shù)值值解解 x x exxy x y yxey )2( 2 1 )(, 1 , 0 1)0( 2 例例6.1 以以 h=0.1為步長，用歐拉法求常微分方程初值問題為步長，用歐拉法求常微分方程初值問題 )9 ,2 , 1 , 0( 1 )( 0 1 i y yexhyy i x iii i 解解：由由歐歐拉拉公公式式得得 xiyiy(xi)y(xi)-yi 0110 0.10.9000000.9093629.362x10-3 . 立表立表 (見表見表6-1（書（書121頁）頁）) 6.1.2 梯形公式與改進(jìn)歐拉公式梯形公式與改進(jìn)歐拉公式

6、 3.梯形公式梯形公式 )56( )( ),(),( 2 )(,()(,( 2 )()( )(,()(,( 2 )(,( 00 111 111 11 1 xyy yxfyxf h yy xyxfxyxf h xyxy xyxfxyxf h dxxyxf iiiiii iiiiii iiii x x i i 從而導(dǎo)出梯形公式從而導(dǎo)出梯形公式 得得 由由 -梯形公式梯形公式也是隱式單步法公式也是隱式單步法公式 ox y ),(yxfY 0 x 1 x 1 i xbxn i x h 圖圖1 梯形公式梯形公式 用梯形公式計(jì)算時(shí)，通常取歐拉公式的解作為迭代初值進(jìn)行迭用梯形公式計(jì)算時(shí)，通常取歐拉公式的解作

7、為迭代初值進(jìn)行迭 代計(jì)算，即采用下式代計(jì)算，即采用下式 ),1 ,0( ),(),( 2 ),( )( 11 )1( 1 )0( 1 k yxfyxf h yy yxfhyy k iiiii k i iiii )1( 1 )( 1 )2( )3( 1 )2( )2( 1 )2( )1( 1 )2( )0( 1 )1( 0 kk yyyyyyy (1) (2) )1( 2 )( 2 )2( )3( 2 )2( )2( 2 )2( )1( 2 )2( )0( 2 )1( )0( 1 kk yyyyyyy 11 )( 1 )1( 1 )( 1 )( ,|, ii k i k i k i yxy yy

8、yk 的的近近似似值值值值 作作為為函函數(shù)數(shù)取取時(shí)時(shí)使使得得迭迭代代計(jì)計(jì)算算到到某某個(gè)個(gè) ) 1, 1 , 0( )2() ,(),( 2 ) 1 (),( , 111 1 1 )1( 1 Ni yxfyxf h yy yxfhyy yy iiiiii iiii ii 則則得得到到顯顯式式計(jì)計(jì)算算公公式式代代替替取取直直接接簡簡單單地地 4.改進(jìn)歐拉公式改進(jìn)歐拉公式 這稱為這稱為改進(jìn)歐拉公式改進(jìn)歐拉公式 NN yyyyyyyyy )2( 3 )2( 3 )1( 2 )2( 2 )1( 1 )2( 1 )1( 0 例例6.2 仍取步長仍取步長h = 0.1，采用改進(jìn)歐拉法重新計(jì)算例，采用改進(jìn)歐拉

9、法重新計(jì)算例 6.1 的的 常微分方程初值問題。常微分方程初值問題。 )9 ,1 , 0( 1 ) ()( 2 )( 0 111 1 1 i y yexyex h yy yexhyy i x ii x iii i x iii ii i (見表見表6-2（書（書125頁）頁）) 這時(shí)改進(jìn)歐拉公式為這時(shí)改進(jìn)歐拉公式為解解 xiyiy(xi)y(xi)-yi 0110 0.10.9095240.9093621.363x10-4 . 立表立表 6.2 計(jì)算公式的誤差分析計(jì)算公式的誤差分析 定義定義6.1 若若 yi+1 是是 yi=y(xi) 從計(jì)算得到的近似解，則稱從計(jì)算得到的近似解，則稱 y(xi

10、+1) yi+1為所用公式的為所用公式的局部截?cái)嗾`差局部截?cái)嗾`差，簡稱為，簡稱為截?cái)嗾`差截?cái)嗾`差。 定理定理6.1 若單步法 yi+1 = yi+h (xi , yi , h) 的局部截?cái)嗾` 差為 O (h p+1) ，且增量函數(shù) (x , y , h) 關(guān)于 y 滿足李普希茲 條件，即存在常數(shù) L0，使對(duì) 成立不等式y(tǒng)y , | | ), ,(),(|yyLhyxhyx 則其整體截?cái)嗾`差 y(xi) yi=O(hp) 截?cái)嗾`差的估計(jì)截?cái)嗾`差的估計(jì)(基本假設(shè)：基本假設(shè)： yi = y( xi ) ) 設(shè)設(shè) y(x) C 3 x0 , b , 則則 )76()()( 2 )()()()( 3 2

11、 1 hOxy h xyhxyhxyxy iiiii （1）對(duì)歐拉公式，有）對(duì)歐拉公式，有 )86()()()( 2 )( )()()(,()( 23 2 11 1 hOhOxy h yxy xyhxyxyxfhxyy iii iiiiii 故 因此，歐拉公式的局部截?cái)嗾`差為 O (h2) （2）對(duì)后退歐拉公式，有）對(duì)后退歐拉公式，有 )96()()()( 2 )( )(,()( 23 2 11 111 hOhOxy h yxy xyxfhxyy iii iiii 故 因此，后退歐拉公式的局部截?cái)嗾`差為 O (h2) （3）對(duì)梯形公式，注意到其公式可改寫為）對(duì)梯形公式，注意到其公式可改寫為 )

12、,(),( 2 1 111 iiiiiii yxfhyyxhfyy 故由式（故由式（6-8）和（）和（6-9）得）得 )106()( )()( 2 )()( 2 2 1 ),()( ),()( 2 1 )( 3 3 2 3 2 111 111 hO hOxy h hOxy h yxfhyxy yxfhyxyyxy ii iiii iiiiii 因此，梯形公式的局部截?cái)嗾`差為因此，梯形公式的局部截?cái)嗾`差為 O ( h3 ) （4）對(duì)改進(jìn)歐拉公式，有）對(duì)改進(jìn)歐拉公式，有 ) ,(),( 2 )( )()( 111 1 iiiiii iii yxfyxf h xyy xyhxyy 而由 ，故有 ),

13、(),(),(yxfyyxfyyxfy yx 得 )()( 2 )()( )()()( )()(,( )()(,() ,( 3 2 1 2 2 )(,( 11 hOxy h xyhxyy hOxyhxy hOfyhfhxyxf xyhxyhxfyxf iiii ii xyxyxii iiiii ii 所以所以 與式（與式（6-7）比較得）比較得 y(xi+1) yi+1 = O ( h3 ) 因此，改進(jìn)歐拉公式的局部截?cái)嗾`差為因此，改進(jìn)歐拉公式的局部截?cái)嗾`差為 O ( h3 ) 定義定義6.2 若一種求解常微分方程初值問題的數(shù)值計(jì)算方法若一種求解常微分方程初值問題的數(shù)值計(jì)算方法 的局部截?cái)嗾`差

14、為的局部截?cái)嗾`差為 O ( hp+1 ) ，則稱該方法為，則稱該方法為 p階精度階精度，或稱該，或稱該 方法為方法為 p階方法階方法。 由此定義知，歐拉方法與后退歐拉方法為一階精度，梯由此定義知，歐拉方法與后退歐拉方法為一階精度，梯 形法與改進(jìn)歐拉方法為二階精度。形法與改進(jìn)歐拉方法為二階精度。 6.3 龍格龍格-庫塔方法庫塔方法 由中值定理，有由中值定理，有 ),(,)(,( )()()()( 1 11 ii iiii xxyfh yxxxyxy 因此，以上介紹的各種單步法本質(zhì)上都是對(duì)平均斜因此，以上介紹的各種單步法本質(zhì)上都是對(duì)平均斜 率率 f( , y( ) 進(jìn)行近似，龍格進(jìn)行近似，龍格-庫

15、塔據(jù)之提出了適當(dāng)選取若庫塔據(jù)之提出了適當(dāng)選取若 干點(diǎn)上的斜率值作近似以構(gòu)造高精度計(jì)算公式的方法，其干點(diǎn)上的斜率值作近似以構(gòu)造高精度計(jì)算公式的方法，其 基本思想是基于泰勒展式的待定系數(shù)法。基本思想是基于泰勒展式的待定系數(shù)法。 6.3.1 二二階階Rung-KuttaRung-Kutta公式公式 問題：問題：建立二階精度的計(jì)算格式形為建立二階精度的計(jì)算格式形為 待定系數(shù)待定系數(shù)為為這里這里ba bhKyahxfK yxfK KKhyy ii ii ii , )126( ),( ),( )( 21 12 1 22111 在在 y(xi) = yi 的假設(shè)下，有的假設(shè)下，有 )()( )()(,( )

16、(,( )()(,( 2 )(,( 2 )(,(1 12 1 hOfybfahxy hOfbhKfahxyxf bhKxyahxfK xyxyxfK ii ii xyxyxi xyxyxii ii iii )()()()( 3 )(,(22 2 211 hOfybfahxyhxyy ii xyxyxiii 故故 解解 )()()()( 3 )(,(22 2 211 hOfybfahxyhxyy ii xyxyxiii 變形歐拉公式變形歐拉公式 根據(jù)格式為二階精度，即根據(jù)格式為二階精度，即 y(xi+1) yi+1 = O(h3) 比較兩式系數(shù)得比較兩式系數(shù)得 )136( 2/1 2/1 1 2

17、 2 21 b a 而而 )( 2 1 2 1 )()( )()()2/()()()( 3 )(,( 2 32 1 hOfyfhxyhxy hOxyhxyhxyxy ii xyxyxii iiii 系數(shù)滿足系數(shù)滿足(6-13)的形為的形為(6-12)計(jì)算格式統(tǒng)稱為二階計(jì)算格式統(tǒng)稱為二階R-K公式。公式。 當(dāng)令當(dāng)令 1=1/2時(shí)，解得 時(shí)，解得 2=1/2 ， ，a=b=1，即為改進(jìn)歐拉公式。，即為改進(jìn)歐拉公式。 若令若令 1=0，解得 ，解得 2=1， ，a=b=1/2，則得另一計(jì)算公式，則得另一計(jì)算公式 )146( )2/, 2/( ),( 12 1 21 hKyhxfK yxfK hKyy

18、 ii ii ii 6.3.2 四階四階 R-K R-K 公式公式 每一步需計(jì)算的 f 值的個(gè)數(shù) 1234567n 8 精度階1234456n-2 1965年，年，Butcher研究發(fā)現(xiàn)顯式研究發(fā)現(xiàn)顯式R-K公式的精度與需要組公式的精度與需要組 合的斜率值的個(gè)數(shù)具有如下關(guān)系合的斜率值的個(gè)數(shù)具有如下關(guān)系 可見，超過四階精度的可見，超過四階精度的R-K公式效率并不高，實(shí)際計(jì)算通公式效率并不高，實(shí)際計(jì)算通 常選用如下四階格式常選用如下四階格式 )156( ),( )2/, 2/( )2/, 2/( ),( )22( 6 34 23 12 1 43211 hKyhxfK hKyhxfK hKyhxfK

19、 yxfK KKKK h yy ii ii ii ii ii 經(jīng)典經(jīng)典R-KR-K公式公式 這時(shí)經(jīng)典這時(shí)經(jīng)典R-K公式為公式為 例例6.3 取步長取步長h = 0.2，采用經(jīng)典，采用經(jīng)典R-K法計(jì)算例法計(jì)算例 6.1 的常微的常微 分方程初值問題。分方程初值問題。 )()( ) 2 () 2 ( ) 2 () 2 ( )22( 6 3 )( 4 2 )2/( 3 1 )2/( 2 1 43211 hKyehxK K h ye h xK K h ye h xK yexK KKKK h yy i hx i i hx i i hx i i x i ii i i i i 取取 h=0.2 計(jì)算得到表計(jì)算得到表6-4(書書133頁）。頁）。 與例與例6.1和例和例6.2比較可見，用經(jīng)典比較可見，用經(jīng)典R-K法計(jì)算得到的解比法計(jì)算得到的解比 用歐拉法和用歐拉法和改進(jìn)歐拉法改進(jìn)歐拉法所得到的解精