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

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

2、xyxy N xb hniihxx bx 6.1 歐拉方法歐拉方法 6.1.1 歐拉公式與改進歐拉公式歐拉公式與改進歐拉公式 算法：算法： N xb hniihxx bx i 0 0 0 ), 1 ,0( , 生生成成節(jié)節(jié)點點作作等等距距分分割割對對區(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 即即 ：積積分分上上對對在在 )46()(,()()( 1

3、1 i i x x ii dxxyxfxyxy 化化為為積積分分方方程程問問題題從從而而把把微微分分方方程程問問題題轉轉 選擇不同的數(shù)值積分公式來求選擇不同的數(shù)值積分公式來求 近似值就得到初值問題的各種數(shù)值解法近似值就得到初值問題的各種數(shù)值解法 1 )(,( i i x x dxxyxf 1.歐拉公式歐拉公式 )(,()()( 1iiii xyxfhxyxy 得得，時時由由 當當 )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 由由此此可可建建立立計計算算公公式式 這稱為這稱為歐拉公式歐拉公式 2.后退歐拉公式后退歐拉公式 ox y ),(yxfY 0 x 1 x 1 i xbxn i x h 得得，時時由由 當當 )46( )(,()(,( 11 1 ii x x xyxfhdxxyxf i i )36()1, 1 , 0( )( ),( 00 111 Ni xyy yxfhyy iiii 由由此此可可建建立立另另一一公公式式 )(,()()( 111 iiii xyxfhxyxy 這稱為這稱為后退歐拉公式后退歐拉公式 后退歐拉公式是一個隱式公式，通常采用迭代法求解。后退歐拉公式是一個隱式公

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 梯形公式與改進歐拉公式梯形公式與改進歐拉公式

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 從而導出梯形公式從而導出梯形公式 得得 由由 -梯形公式梯形公式也是隱式單步法公式也是隱式單步法公式 ox y ),(yxfY 0 x 1 x 1 i xbxn i x h 圖圖1 梯形公式梯形公式 用梯形公式計算時，通常取歐拉公式的解作為迭代初值進行迭用梯形公式計算時，通常取歐拉公式的解作

7、為迭代初值進行迭 代計算，即采用下式代計算，即采用下式 ),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ù)取取時時使使得得迭迭代代計計算算到到某某個個 ) 1, 1 , 0( )2() ,(),( 2 ) 1 (),( , 111 1 1 )1( 1 Ni yxfyxf h yy yxfhyy yy iiiiii iiii ii 則則得得到到顯顯式式計計算算公公式式代代替替取取直直接接簡簡單單地地 4.改進歐拉公式改進歐拉公式 這稱為這稱為改進歐拉公式改進歐拉公式 NN yyyyyyyyy )2( 3 )2( 3 )1( 2 )2( 2 )1( 1 )2( 1 )1( 0 例例6.2 仍取步長仍取步長h = 0.1，采用改進歐拉法重新計算例，采用改進歐拉

9、法重新計算例 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頁）頁）) 這時改進歐拉公式為這時改進歐拉公式為解解 xiyiy(xi)y(xi)-yi 0110 0.10.9095240.9093621.363x10-4 . 立表立表 6.2 計算公式的誤差分析計算公式的誤差分析 定義定義6.1 若若 yi+1 是是 yi=y(xi) 從計算得到的近似解，則稱從計算得到的近似解，則稱 y(xi

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

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

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 因此，梯形公式的局部截斷誤差為因此，梯形公式的局部截斷誤差為 O ( h3 ) （4）對改進歐拉公式，有）對改進歐拉公式，有 ) ,(),( 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 ) 因此，改進歐拉公式的局部截斷誤差為因此，改進歐拉公式的局部截斷誤差為 O ( h3 ) 定義定義6.2 若一種求解常微分方程初值問題的數(shù)值計算方法若一種求解常微分方程初值問題的數(shù)值計算方法 的局部截斷誤差

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

15、塔據(jù)之提出了適當選取若庫塔據(jù)之提出了適當選取若 干點上的斜率值作近似以構造高精度計算公式的方法，其干點上的斜率值作近似以構造高精度計算公式的方法，其 基本思想是基于泰勒展式的待定系數(shù)法?；舅枷胧腔谔├照故降拇ㄏ禂?shù)法。 6.3.1 二二階階Rung-KuttaRung-Kutta公式公式 問題：問題：建立二階精度的計算格式形為建立二階精度的計算格式形為 待定系數(shù)待定系數(shù)為為這里這里ba bhKyahxfK yxfK KKhyy ii ii ii , )126( ),( ),( )( 21 12 1 22111 在在 y(xi) = yi 的假設下，有的假設下，有 )()( )()(,( )

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)計算格式統(tǒng)稱為二階計算格式統(tǒng)稱為二階R-K公式。公式。 當令當令 1=1/2時，解得 時，解得 2=1/2 ， ，a=b=1，即為改進歐拉公式。，即為改進歐拉公式。 若令若令 1=0，解得 ，解得 2=1， ，a=b=1/2，則得另一計算公式，則得另一計算公式 )146( )2/, 2/( ),( 12 1 21 hKyhxfK yxfK hKyy

18、 ii ii ii 6.3.2 四階四階 R-K R-K 公式公式 每一步需計算的 f 值的個數(shù) 1234567n 8 精度階1234456n-2 1965年，年，Butcher研究發(fā)現(xiàn)顯式研究發(fā)現(xiàn)顯式R-K公式的精度與需要組公式的精度與需要組 合的斜率值的個數(shù)具有如下關系合的斜率值的個數(shù)具有如下關系 可見，超過四階精度的可見，超過四階精度的R-K公式效率并不高，實際計算通公式效率并不高，實際計算通 常選用如下四階格式常選用如下四階格式 )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公式公式 這時經(jīng)典這時經(jīng)典R-K公式為公式為 例例6.3 取步長取步長h = 0.2，采用經(jīng)典，采用經(jīng)典R-K法計算例法計算例 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 計算得到表計算得到表6-4(書書133頁）。頁）。 與例與例6.1和例和例6.2比較可見，用經(jīng)典比較可見，用經(jīng)典R-K法計算得到的解比法計算得到的解比 用歐拉法和用歐拉法和改進歐拉法改進歐拉法所得到的解精