常微分方程數(shù)值解法

1、數(shù)數(shù) 值值 分分 析析 第5章 常微分方程數(shù)值解法 1 引言 v1.0 基本概念基本概念 1. 常微分方程的初值問題： 稱為具有初值稱為具有初值(1.2)的常微分方程的常微分方程. 若若f(x,y)在在a x b, |y|+ 上連續(xù)，且關(guān)于上連續(xù)，且關(guān)于y滿足滿足Lip 條件：條件： 常數(shù)常數(shù)L使使| f(x, y1) f(x, y2)| L|y1 y2| 則初值問題則初值問題(1.1)(1.2)存在唯一連續(xù)可微解存在唯一連續(xù)可微解y(x). 注：以下總假設(shè)注：以下總假設(shè)f 滿足滿足Lip條件條件. )2 . 1()( )1 . 1(,),( 0 yay baxyxfy 1 引言 v1.0 基

2、本概念基本概念 1. 常微分方程的初值問題： 稱為具有初值稱為具有初值(1.2)的常微分方程的常微分方程. (1.1)(1.2)等價(jià)于微分方程：等價(jià)于微分方程： (1.3) 注：一般無初等解注：一般無初等解(解析解解析解)，即使有形式也復(fù)雜，即使有形式也復(fù)雜. )2 . 1()( )1 . 1(,),( 0 yay baxyxfy .)(,()( 0 x a dttytfyxy 1 引言 v1.0 基本概念基本概念 2. 初值問題的數(shù)值解 設(shè)設(shè)(1.1)(1.2)的解的解y(x)在節(jié)點(diǎn)在節(jié)點(diǎn)xi處的近似解值為處的近似解值為 yi y(xi), a x1 x2 xn = b 則稱則稱yi (i

3、= 1, 2, , n)為為(1.1)(1.2)的數(shù)值解，又稱的數(shù)值解，又稱y(xi) 的計(jì)算值的計(jì)算值. )2 . 1()( )1 . 1(,),( 0 yay baxyxfy 1 引言 v1.0 基本概念基本概念 3. 數(shù)值方法 兩種轉(zhuǎn)化：兩種轉(zhuǎn)化： 由微分出發(fā)的數(shù)值方法由微分出發(fā)的數(shù)值方法. 由積分由積分 出發(fā)的數(shù)值方法出發(fā)的數(shù)值方法. 計(jì)算方法計(jì)算方法 步進(jìn)法：從初始條件出發(fā)，逐步求步進(jìn)法：從初始條件出發(fā)，逐步求y1, y2, , yn. 又有兩種：單步法，多步法又有兩種：單步法，多步法. 注：采用等距節(jié)點(diǎn)：注：采用等距節(jié)點(diǎn)： x a dttytf)(,( . ),.,2 , 1 ,

4、0(.ni n ab hihaxi 1 引言 v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式. (1.6) ).( 2 )()( 1 )( ).( 2 )()( 1 )( 11 1 iiii iiii y h xyxy h xy y h xyxy h xy ,)( 2 )( )( 2 )()()( 1 11 1 iiiiii iiii xxy h xy y h xyxyxy h 1 引言 v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式. 1. 前進(jìn)歐拉公式 (1.6)的前半部分為：的前半部分為： 令令 yi+1 = yi + hf(xi, yi) (1.7) 其中其中yi =

5、y(xi) , 則則yi+1 y(xi+1) ).( 2 )(,()()( 1 1iiiii y h xyxfxyxy h ).( 2 )(,()()( 2 1iiiii y h xyxhfxyxy ,)( 2 )( )( 2 )()()( 1 11 1 iiiiii iiii xxy h xy y h xyxyxy h 1 引言 v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式. 1. 前進(jìn)歐拉公式 令令 yi+1 = yi + hf(xi, yi) (1.7) 其中其中yi = y(xi) , 則則yi+1 y(xi+1) 記記 (1.8) 則則 稱稱(1.7)為前進(jìn)歐拉求解公式為前

6、進(jìn)歐拉求解公式. 簡稱為歐拉公式或歐拉法簡稱為歐拉公式或歐拉法. (1.8)稱為歐拉公式的余項(xiàng)：稱為歐拉公式的余項(xiàng)：ei+1(h) = y(xi+1) yi+1 ).( 2 )(,()()( 2 1iiiii y h xyxhfxyxy )( 2 )( 2 1ii y h he ).( 2 )( 2 1ii xy h he 1 引言 v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式. 2. 后退歐拉公式 (1.6)的后半部分的后半部分 令令 yi+1 = yi + hf(xi+1, yi+1) (1.9) 其中其中yi = y(xi), 則則yi+1 y(xi+1) ).( 2 )()(

7、)( 1 11iiii y h xyxyxy h ).( 2 )(,()()( 2 111iiiii y h xyxhfxyxy ,)( 2 )( )( 2 )()()( 1 11 1 iiiiii iiii xxy h xy y h xyxyxy h 1 引言 v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式. 2. 后退歐拉公式 令令 yi+1 = yi + hf(xi+1, yi+1) (1.9) 其中其中yi = y(xi), 則則yi+1 y(xi+1) 注：注：(1.9)中中f(xi+1, yi+1) f(xi+1, y(xi+1) 余項(xiàng)余項(xiàng) (1.10) ).( 2 )(

8、2 ),()(,()()( 1111111 i iiiiiiii y h y h yxfxyxfhyxyhe 1 引言 v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式. 2. 后退歐拉公式 令令 yi+1 = yi + hf(xi+1, yi+1) (1.9) 其中其中yi = y(xi), 則則yi+1 y(xi+1) 注：注： 稱稱(1.9)為后退歐拉公式為后退歐拉公式(后退歐拉法后退歐拉法). 稱稱(1.10)為后退歐拉法的誤差近似值為后退歐拉法的誤差近似值. 歐拉法與后退歐拉公式的區(qū)別：歐拉法與后退歐拉公式的區(qū)別： (1.7)為直接計(jì)算公式稱顯式公式為直接計(jì)算公式稱顯式公式.

9、(1.9)為關(guān)于函數(shù)方程稱隱式公式為關(guān)于函數(shù)方程稱隱式公式. 1 引言 v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式. 【例例1】取取h=0.1求解初值問題：求解初值問題： (1.11). 解：解： ，xi = ih = 0.1 i, (i = 0, 1, 2, , 10) 歐拉法：歐拉法： 1)0( 1 , 0 2 y x y x yy y x yyxf 2 ),( . 9,.,2 , 1 , 0) 2 ( 1 i y x yhyy i i iii 1 引言 v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式. 【例例1】取取h=0.1求解初值問題：求解初值問題： (1.11)

10、. 解：解： ，xi = ih = 0.1 i, (i = 0, 1, 2, , 10) 后退歐拉法：后退歐拉法： 1)0( 1 , 0 2 y x y x yy y x yyxf 2 ),( ) 2 ( 1 1 11 i i iii y x yhyy . 9,.,2 , 1 , 0 )1(2 )1(8 1 2 1 i h hxhyy y iii i 1 引言 v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式. 注：為避免求解函數(shù)方程，采用顯式與隱式結(jié)合的方注：為避免求解函數(shù)方程，采用顯式與隱式結(jié)合的方 法：法： 此方法稱為此方法稱為 預(yù)測預(yù)測校正系統(tǒng)校正系統(tǒng). 求解過程為：求解過程為：

11、 校校正正值值隱隱式式 預(yù)預(yù)測測值值顯顯式式 )( )( ),( ),( 111 1 iiii iiii yxhfyy yxhfyy .)()()( 22110nn yyyyyyy 1 引言 v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式. 預(yù)測預(yù)測校正系統(tǒng)：校正系統(tǒng)： 【例例2】利用預(yù)測利用預(yù)測校正系統(tǒng)求解例校正系統(tǒng)求解例1. 校校正正值值隱隱式式 預(yù)預(yù)測測值值顯顯式式 )( )( ),( ),( 111 1 iiii iiii yxhfyy yxhfyy 1)0( 1 , 0 2 y x y x yy . 9,.,2 , 1 , 0 ) 2 ( ) 2 ( 1 1 11 1 i y

12、 x yhyy y x yhyy i i iii i i iii 1 引言 v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式. 預(yù)測預(yù)測校正系統(tǒng)：校正系統(tǒng)： 注：顯式比隱式方便，但有時(shí)隱式效果比顯式好注：顯式比隱式方便，但有時(shí)隱式效果比顯式好.(4 介紹介紹). 校校正正值值隱隱式式 預(yù)預(yù)測測值值顯顯式式 )( )( ),( ),( 111 1 iiii iiii yxhfyy yxhfyy 1 引言 v1.2 截?cái)嗾`差截?cái)嗾`差 定義定義1.1 稱稱ek(h) = y(xk) yk為計(jì)算為計(jì)算yk的公式第的公式第k步的局步的局 部截?cái)嗾`差部截?cái)嗾`差. 注：注：“局部局部”是指在計(jì)算第是指

13、在計(jì)算第k步時(shí)，假定前面步時(shí)，假定前面yi = y(xi) (i k).而而yk y(xk) 歐拉法歐拉法. 后退歐拉法后退歐拉法. 一般根據(jù)一般根據(jù)y(xk)對(duì)對(duì)y( k), y( k)做估計(jì)做估計(jì). )( 2 )( 2 1kk y h he ).( 2 )( 2 1ii y h he 1 引言 v1.2 截?cái)嗾`差截?cái)嗾`差 定義定義1.2 設(shè)設(shè)ei(h) (i = 1, 2, , k)為求解公式第為求解公式第i步的局部步的局部 截?cái)嗾`差截?cái)嗾`差.稱稱 為該求解公式在點(diǎn)上的整體截?cái)嗾`差為該求解公式在點(diǎn)上的整體截?cái)嗾`差. 注：局部截?cái)嗾`差注：局部截?cái)嗾`差ek(h)與與yk有關(guān)有關(guān). 整體截?cái)嗾`差

14、整體截?cái)嗾`差Ek(h)與與y1, y2, , yk有關(guān)有關(guān). 所有所有ek(h)都與都與h有關(guān)有關(guān). k i ik hehE 1 )()( 1 引言 v1.2 截?cái)嗾`差截?cái)嗾`差 定義定義1.3 若局部截?cái)嗾`差若局部截?cái)嗾`差e(h)=O(hp+1)，則稱該求解公，則稱該求解公 式具有式具有p階精度階精度. 注：歐拉法具有一階精度注：歐拉法具有一階精度.(精度越高越好精度越高越好) 1 引言 作業(yè)作業(yè) P208 1，2，3. 1 引言 v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 (1.13) 若已知若已知y(xk) = yk, 則計(jì)算積分可求出則計(jì)算積分可求出y(xk+1) . 如用矩

15、形公式求積分如用矩形公式求積分 則有則有y(xk+1) = y(xk) + hf(xk, yk) 令令yk+1 = y(xk) + hf(xk, yk)即為歐拉公式即為歐拉公式. 故歐拉公式又故歐拉公式又 稱矩形法稱矩形法. ).()()()(,( 1 11 kk x x x x xyxydxxydxxyxf k k k k .)(,()()( 1 1 k k x x kk dxxyxfxyxy ).,()(,( 1 kk x x xxhfdxxyxf k k 1 引言 v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 (1.13) 考慮考慮 1. 梯形公式 記記 (1.14) .)(,

16、()()( 1 1 k k x x kk dxxyxfxyxy 11 )()(,( k k k k x x x x dxxFdxxyxf ).(,()(,( 2 )( )(,()(,( 2 )(,( 111 11 1 kkkkkk kkkk x x xyxfxyxf h yxy xyxfxyxf h dxxyxf k k ),(),( 2 111 kkkkkk yxfyxf h yy 1 引言 v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 1. 梯形公式 記記 (1.14) 稱稱(1.14)為梯形為梯形(求解求解)公式公式. 簡稱梯形法簡稱梯形法. ).(,()(,( 2 )( )(

17、,()(,( 2 )(,( 111 11 1 kkkkkk kkkk x x xyxfxyxf h yxy xyxfxyxf h dxxyxf k k ),(),( 2 111 kkkkkk yxfyxf h yy 1 引言 v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 1. 梯形公式 梯形梯形(求解求解)公式公式, 簡稱梯形法簡稱梯形法: (1.14) 注：梯形公式的余項(xiàng)：注：梯形公式的余項(xiàng)： 故是二階精度故是二階精度. ),(),( 2 111 kkkkkk yxfyxf h yy ,)( 12 )( )( 1 33 1 iii i i xxhOh y he 3 )( 12 )(

18、 ab f RT v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 1. 梯形公式 (1.14) 梯形公式為隱式公式梯形公式為隱式公式. 預(yù)測預(yù)測校正系統(tǒng)校正系統(tǒng) (1.15) 稱稱(1.15)為改進(jìn)的歐拉公式，也可記為為改進(jìn)的歐拉公式，也可記為 1 引言 ),(),( 2 111 kkkkkk yxfyxf h yy . 1,.,2 , 1 , 0 )( )( ),( ),( 111 1 ni yxhfyy yxhfyy iiii iiii 梯形梯形校正校正 歐拉歐拉預(yù)測預(yù)測 ),(,(),( 2 11iiiiiiii yxhfyxfyxf h yy 1 引言 v1.3 基于數(shù)值積分的

19、求解公式基于數(shù)值積分的求解公式 1. 梯形公式 (1.14) 可以證明，改進(jìn)歐拉公式也具有二階精度可以證明，改進(jìn)歐拉公式也具有二階精度. ),(),( 2 111 kkkkkk yxfyxf h yy 1 引言 v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 【例例3】用歐拉法，梯形法以及改進(jìn)歐拉法求解用歐拉法，梯形法以及改進(jìn)歐拉法求解 取取h=0.1.計(jì)算到計(jì)算到x=0.5. 解：解：f(x, y) = xy + 1, a = x0 = 0, b = 0.5, y0 = 1, n = 5 (Euler法法) 求解公式：求解公式：yk =yk1+h(xk1yk1+1) = hxk1+(

20、1 h)yk1 + h = 0.1xk1+0.9yk1+0.1 1)0( 1 y xyy 1 引言 v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 【例例3】用歐拉法，梯形法以及改進(jìn)歐拉法求解用歐拉法，梯形法以及改進(jìn)歐拉法求解 解：解：f(x, y) = xy + 1, a = x0 = 0, b = 0.5, y0 = 1, n = 5 (梯形法梯形法)求解公式：求解公式： yk=yk1+h(xk1yk1+1)+(xkyk+1)/2 解出解出yk，得，得 1)0( 1 y xyy 1 . 2 2 . 0)( 1 . 09 . 1 2 2)()2( 11 11 kkk kkk k xx

21、y h hxxhyh y 方程方程 1 引言 v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 【例例3】用歐拉法，梯形法以及改進(jìn)歐拉法求解用歐拉法，梯形法以及改進(jìn)歐拉法求解 解：解：f(x, y) = xy + 1, a = x0 = 0, b = 0.5, y0 = 1, n = 5 (改進(jìn)改進(jìn)Euler法法)求解公式：求解公式： yk=yk1+h(xk1yk1+1) + xk (yk +h(xkyk+1)+1/2 得得 =0.905yk1+0.045xk1+0.05xk+0.095 1)0( 1 y xyy 2 )2( 22 )1( ) 2 )2( 1( 11 hh x h x h

22、h y hh kkk 方程方程 1 引言 v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 2. 辛卜生公式 記記 (1.17) )(,() 2 (, 2 (4)(,( 6 )(,( 11 11 1 kk kkkk kk x x xyxf xx y xx fxyxf h dxxyxf k k )(,() 2 (, 2 (4)(,( 6 1111 hxyhxf h xy h xfxyxf h kkkkkk )(,( ) 2 (, 2 (4)(,( 6 11 111 hxyhxf h xy h xfxyxf h yy kk kkkkkk 1 引言 v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分

23、的求解公式 2. 辛卜生公式 記記 (1.17) 其余項(xiàng)其余項(xiàng) )(,( ) 2 (, 2 (4)(,( 6 11 111 hxyhxf h xy h xfxyxf h yy kk kkkkkk )(,( 16180 )( )4( 5 kkk yf h he ,)()( 16180 1 5)5( 5 iiii xxhOy h 1 引言 v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 2. 辛卜生公式 記記 (1.17) 將將xk1, xk 對(duì)分：對(duì)分： 調(diào)整下標(biāo)為調(diào)整下標(biāo)為xi2, xi ： xi2 = xk1, xi1 = xk1+h1, xi = xk1+2h1= xk 則則(1.

24、17)化為化為 (1.19) 稱稱(1.19)為辛卜生求解公式，其中為辛卜生求解公式，其中fk2= f(xk2, y(xk2)， fk1 = f(xk1, y(xk1)，fk = f(xk, y(xk) )(,( ) 2 (, 2 (4)(,( 6 11 111 hxyhxf h xy h xfxyxf h yy kk kkkkkk . 2 1 h h 4 3 12 1 2iiiii fff h yy 1 引言 v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 2. 辛卜生公式 記記 (1.17) (1.19) 稱稱(1.19)為辛卜生求解公式，其中為辛卜生求解公式，其中fi2= f(x

25、i2, y(xi2)， fi1 = f(xi1, y(xi1)，fi = f(xi, y(xi) 注：注： (1.19)的誤差：的誤差： )(,( ) 2 (, 2 (4)(,( 6 11 111 hxyhxf h xy h xfxyxf h yy kk kkkkkk 4 3 12 1 2iiiii fff h yy );()( 90 )( 16180 )2( )( 5 1 )( 5 1)( 5 1 1 hOy h y h he k k k k k 1 引言 v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 2. 辛卜生公式 記記 (1.17) (1.19) 稱稱(1.19)為辛卜生求解

26、公式，其中為辛卜生求解公式，其中fi2= f(xi2, y(xi2)， fi1 = f(xi1, y(xi1)，fi = f(xi, y(xi) 注：注： 隱式隱式(需顯化需顯化) 多步多步將在將在3中討論中討論. )(,( ) 2 (, 2 (4)(,( 6 11 111 hxyhxf h xy h xfxyxf h yy kk kkkkkk 4 3 12 1 2iiiii fff h yy 2 Runge - Kutta法 v2.0 原理原理 其中其中K = f( , y( ) = y( )稱為稱為y在在xi1, xi上的平均斜率上的平均斜率. 歐拉法：歐拉法： 改進(jìn)歐拉法：改進(jìn)歐拉法：

27、(2.1) )(,()()(,()()( 11 1 yhfxydxxyxfxyxy i x x ii i i hKxy i )( 1 ),( 111 11 ii ii yxfK hKyy ),( ),( ) 2 ( 112 111 21 1 hKyxfK yxfK KK hyy ii ii ii 2 Runge - Kutta法 v2.0 原理原理 其中其中K = f( , y( ) = y( )稱為稱為y在在xi1, xi上的平均斜率上的平均斜率. 對(duì)對(duì)(1.17)顯化：顯化： 辛卜生：辛卜生： (2.4) )(,()()(,()()( 11 1 yhfxydxxyxfxyxy i x x

28、ii i i hKxy i )( 1 )2(,( ) 2 , 2 ( ),( 6 4 12113 1112 111 321 1 KKhyhxfK K h y h xfK yxfK KKK hyy ii ii ii ii 2 Runge - Kutta法 v2.0 原理原理 其中其中K = f( , y( ) = y( )稱為稱為y在在xi1, xi上的平均斜率上的平均斜率. 設(shè)想：在中多計(jì)算設(shè)想：在中多計(jì)算(預(yù)測預(yù)測)幾個(gè)點(diǎn)上的值然后可加權(quán)取幾個(gè)點(diǎn)上的值然后可加權(quán)取 平均值作為的近似值可能構(gòu)成更高階的公式平均值作為的近似值可能構(gòu)成更高階的公式. hKxyxy ii )()( 1 ),( 111

29、 11 ii ii yxfK hKyy ),( ),( ) 2 ( 112 111 21 1 hKyxfK yxfK KK hyy ii ii ii )2(,( ) 2 , 2 ( ),( 6 4 12113 1112 111 321 1 KKhyhxfK K h y h xfK yxfK KKK hyy ii ii ii ii 一階一階 二階二階 三階三階 2 Runge - Kutta法 v2.1 Runge - Kutta公式公式 (*) 其中其中0 j 1，yi1 + jh是是y(xi1 + jh) 的預(yù)測值的預(yù)測值. 稱稱(*)為為R-K公式公式 注：注：(2.1)(2.4)分別稱為

30、二階，三階分別稱為二階，三階R-K公式公式. j， j， j為待定系數(shù)為待定系數(shù). 使使(*)的階數(shù)盡量高的階數(shù)盡量高. ),( ),( ),( )( 11 21212 111 22111 hyhxfK hyhxfK yxfK KKKhyy mimim ii ii mmii 2 Runge - Kutta法 v2.1 Runge - Kutta公式公式 參數(shù)的確定，以參數(shù)的確定，以m = 2為例為例. 欲求欲求 1， 2， 2 . ),()(,( ),( 1212121212 111 22111 hKyhxfhxyhxfK yxfK KKhyy iiii ii ii 原則原則: 使使ei(h)

31、 = y(xi) yi的階數(shù)盡可能高的階數(shù)盡可能高 2 Runge - Kutta法 v2.1 Runge - Kutta公式公式 展開展開 展開展開 ),()(,( ),( 1212121212 111 22111 hKyhxfhxyhxfK yxfK KKhyy iiii ii ii 原則原則: 使使ei(h) = y(xi) yi的階數(shù)盡可能高的階數(shù)盡可能高 ).()( 2 )()()()( 3 1 2 111 hOxy h xyhxyhxyxy iiiii )(),(),( 2 11122112 hOyxf y hK x hyxfK iiii )(),(),(),( 2 1111121

32、1 hOyxf y Kyxf x hyxf iiiiii 2 Runge - Kutta法 v2.1 Runge - Kutta公式公式 ),()(,( ),( 1212121212 111 22111 hKyhxfhxyhxfK yxfK KKhyy iiii ii ii 原則原則: 使使ei(h) = y(xi) yi的階數(shù)盡可能高的階數(shù)盡可能高 ).()( 2 )()()()( 3 1 2 111 hOxy h xyhxyhxyxy iiiii )(),(),(),( 2 111112112 hOyxf y Kyxf x hyxfK iiiiii . )()()()( )(),(),(

33、),(),( 3 1 2 221211 3 11111 2 22 1121111 hOxyhxyhxy hOyxf y Kyxf x h yxhfyxhfyy iii iiii iiiiii 2 Runge - Kutta法 v2.1 Runge - Kutta公式公式 欲求截?cái)嗾`差欲求截?cái)嗾`差ei(h) = y(xi) yi關(guān)于關(guān)于h的階數(shù)盡可能高，的階數(shù)盡可能高， 應(yīng)使應(yīng)使 ),()(,( ),( 1212121212 111 22111 hKyhxfhxyhxfK yxfK KKhyy iiii ii ii 無窮多解，從而有許多無窮多解，從而有許多2階階R-K公式公式 ).()( 2 )

34、()()()( 3 1 2 111 hOxy h xyhxyhxyxy iiiii ).()()()( 3 1 2 221211 hOxyhxyhxyy iiii 2 1 1 22 21 2 Runge - Kutta法 v2.1 Runge - Kutta公式公式 應(yīng)使應(yīng)使 注：注： 取取 1= 2= 1/2， 2 = 1，即為改進(jìn)歐拉公式，即為改進(jìn)歐拉公式. ),()(,( ),( 1212121212 111 22111 hKyhxfhxyhxfK yxfK KKhyy iiii ii ii 2 1 1 22 21 ),( ),( ) 2 ( 112 111 21 1 hKyxfK yx

35、fK KK hyy ii ii ii 2 Runge - Kutta法 v2.1 Runge - Kutta公式公式 應(yīng)使應(yīng)使 注：注：取取 1= 0， 2 = 1， 2 = 1/2，即為中點(diǎn)公式，即為中點(diǎn)公式 ),()(,( ),( 1212121212 111 22111 hKyhxfhxyhxfK yxfK KKhyy iiii ii ii 2 1 1 22 21 ) 2 , 2 ( ),( 1112 111 21 K h y h xfK yxfK hKyy ii ii ii 2 Runge - Kutta法 v2.1 Runge - Kutta公式公式 應(yīng)使應(yīng)使 注：二階注：二階R-K

36、公式的截?cái)嗾`差為故為二階方法公式的截?cái)嗾`差為故為二階方法.相仿相仿 可得更高階的可得更高階的R-K公式公式. ),()(,( ),( 1212121212 111 22111 hKyhxfhxyhxfK yxfK KKhyy iiii ii ii 2 1 1 22 21 2 Runge - Kutta法 v2.2 經(jīng)典經(jīng)典R-K公式公式 在在4解解R-K公式中最重要的是經(jīng)典公式中最重要的是經(jīng)典R-K公式公式. (2.6) 注：注： (2.6)為為4階方法階方法. ),( ) 2 , 2 ( ) 2 , 2 ( ),( 22 6 3114 2113 1112 111 43211 hKyhxfK

37、K h y h xfK K h y h xfK yxfK KKKK h yy ii ii ii ii ii 2 Runge - Kutta法 v2.2 經(jīng)典經(jīng)典R-K公式公式 在在4解解R-K公式中最重要的是經(jīng)典公式中最重要的是經(jīng)典R-K公式公式. (2.6) 注：注：R-K法對(duì)法對(duì)4階以上不一定能提高整數(shù)階階以上不一定能提高整數(shù)階. ),( ) 2 , 2 ( ) 2 , 2 ( ),( 22 6 3114 2113 1112 111 43211 hKyhxfK K h y h xfK K h y h xfK yxfK KKKK h yy ii ii ii ii ii 2 Runge - K

38、utta法 v2.2 經(jīng)典經(jīng)典R-K公式公式 【例例4】使用三階，四階使用三階，四階R-K法求解初值問題：法求解初值問題： 的部分計(jì)算值的部分計(jì)算值y1，y2，y3，其中，其中h=0.1. 解解 使用三階使用三階R-K法法 1)0( 5 . 0 , 0 2 y xyy 2 12112113 2 111112 2 1111 321 1 )2()2(,( ) 2 () 2 , 2 ( ),( 6 4 KKhyKKhyhxfK K h yK h y h xfK yyxfK KKK hyy iii iii iii ii 2 Runge - Kutta法 【例例4】使用三階，四階使用三階，四階R-K法求

39、解初值問題：法求解初值問題： 的部分計(jì)算值的部分計(jì)算值y1，y2，y3，其中，其中h=0.1. 解解 使用三階使用三階R-K法法 1)0( 5 . 0 , 0 2 y xyy 2 1213 2 112 2 11 321 1 )2( ) 2 ( 6 4 KKhyK K h yK yK KKK hyy i i i ii 2 Runge - Kutta法 【例例4】使用三階，四階使用三階，四階R-K法求解初值問題：法求解初值問題： 的部分計(jì)算值的部分計(jì)算值y1，y2，y3，其中，其中h=0.1. 解解 使用四階使用四階R-K法法 1)0( 5 . 0 , 0 2 y xyy 2 313114 2 2

40、12113 2 111112 2 1111 43211 )(),( ) 2 () 2 , 2 ( ) 2 () 2 , 2 ( )(),( 22 6 hKyhKyhxfK K h yK h y h xfK K h yK h y h xfK yyxfK KKKK h yy iii iii iii iii ii 2 Runge - Kutta法 【例例4】使用三階，四階使用三階，四階R-K法求解初值問題：法求解初值問題： 的部分計(jì)算值的部分計(jì)算值y1，y2，y3，其中，其中h=0.1. 解解 使用四階使用四階R-K法法 1)0( 5 . 0 , 0 2 y xyy 2 314 2 213 2 11

41、2 2 11 43211 )( ) 2 ( ) 2 ( )( 22 6 hKyK K h yK K h yK yK KKKK h yy i i i i ii 2 Runge - Kutta法 注注 使用使用R-K法要求具備較好的光滑性，否則效果不如法要求具備較好的光滑性，否則效果不如 低階的低階的. 作業(yè)作業(yè)P209 8 9，10. 3 線性多步法 單步法的優(yōu)點(diǎn)：簡單，計(jì)算單步法的優(yōu)點(diǎn)：簡單，計(jì)算yk+1只用只用yk. 缺點(diǎn)缺點(diǎn): 沒有充分利用前面的信息且計(jì)算沒有充分利用前面的信息且計(jì)算y(xk+ h)較困難較困難 回顧回顧Simpson： (1.19) 考慮：考慮： (3.1) 兩種插值求積

42、：兩種插值求積： 將將xk1, xk增加內(nèi)部節(jié)點(diǎn)，改為增加內(nèi)部節(jié)點(diǎn)，改為xk2, xk導(dǎo)出的公式導(dǎo)出的公式 稱為閉型求解公式稱為閉型求解公式. 4 3 12 1 2kkkkk fff h yy 線性多步線性多步 .)(,()( 1 1 k k x x kk dxxyxfyxy 3 線性多步法 考慮：考慮： (3.1) 兩種插值求積：兩種插值求積： 將將xk1, xk增加內(nèi)部節(jié)點(diǎn)，改為增加內(nèi)部節(jié)點(diǎn)，改為xk2, xk導(dǎo)出的公式導(dǎo)出的公式 稱為閉型求解公式稱為閉型求解公式. 在在xk1, xk外增加插值節(jié)點(diǎn)，導(dǎo)出的公式稱為開型外增加插值節(jié)點(diǎn)，導(dǎo)出的公式稱為開型 求解公式求解公式. 開型有顯和隱，閉

43、型也有顯和隱開型有顯和隱，閉型也有顯和隱. .)(,()( 1 1 k k x x kk dxxyxfyxy 3 線性多步法 v3.1 開型求解公式開型求解公式 1. 亞當(dāng)斯顯式求解公式 取節(jié)點(diǎn)取節(jié)點(diǎn)xk3, xk2, xk1，在，在xk3, xk上作上作F(x) = f(x, y(x) 的插值多項(xiàng)式的插值多項(xiàng)式. 3 1 3 2313 21 2 3212 31 1 3121 32 22 )( ! 3 )( )( )( )( )( )( )( )( )( )( )()()( i ikk kkkk kk k kkkk kk k kkkk kk xx F xF xxxx xxxx xF xxxx

44、xxxx xF xxxx xxxx xRxLxF 3 線性多步法 v3.1 開型求解公式開型求解公式 1. 亞當(dāng)斯顯式求解公式 取節(jié)點(diǎn)取節(jié)點(diǎn)xk3, xk2, xk1，在，在xk3, xk上上 記記xki = xk ih, x = xk + th，則，則 3 1 3 2313 21 2 3212 31 1 3121 32 22 )( ! 3 )( )( )( )( )( )( )( )( )( )( )()()( i ikk kkkk kk k kkkk kk k kkkk kk xx F xF xxxx xxxx xF xxxx xxxx xF xxxx xxxx xRxLxF . 2 )2

45、)(1( )3)(1( 2 )3)(2( )()( 32122 kkk F tt FttF tt tLxL 3 線性多步法 v3.1 開型求解公式開型求解公式 1. 亞當(dāng)斯顯式求解公式 取節(jié)點(diǎn)取節(jié)點(diǎn)xk3, xk2, xk1，記，記xki = xk ih, x = xk + th，則，則 代入代入(3.1)得得 . 2 )2)(1( )3)(1( 2 )3)(2( )()( 32122 kkk F tt FttF tt tLxL .51623 12 2 )2)(1( )3)(1( 2 )3)(2( )()()( 3211 3 0 1 2 0 1 1 0 1 1 0 1 2121 1 kkkk

46、kkkk k x x kk FFF h y Fdt tt FdtttFdt tt hy dttLhydxxLyxy k k 3 線性多步法 v3.1 開型求解公式開型求解公式 1. 亞當(dāng)斯顯式求解公式 取節(jié)點(diǎn)取節(jié)點(diǎn)xk3, xk2, xk1，記，記xki = xk ih, x = xk + th，則，則 令令 (3.4) 稱稱(3.4)為亞當(dāng)斯顯式求解公式為亞當(dāng)斯顯式求解公式(線性多步線性多步). .51623 12 )( 3211 kkkkk FFF h yxy ).,(5),(16),(23 12 3322111 kkkkkkkk yxFyxFyxF h yy 3 線性多步法 v3.1 開

47、型求解公式開型求解公式 1. 亞當(dāng)斯顯式求解公式 取節(jié)點(diǎn)取節(jié)點(diǎn)xk3, xk2, xk1，記，記xki = xk ih, x = xk + th，則，則 余項(xiàng)：余項(xiàng)： 從而從而(3.4)具有具有3階精度階精度. 稱為稱為3階亞當(dāng)斯求解公式階亞當(dāng)斯求解公式. .)1( ! 3 )( )( ! 3 )( )( 3 3 1 )4(3 1 3 ht y xx F xR ii ik .)( 8 3 )3)(2)(1()( 6 )1()( 6 )( 4)4( 0 1 1 )4( 4 0 1 3 1 )4( 4 3 1 hydtttty h dtty h dxxR i x x k k 3 線性多步法 v3.

48、1 開型求解公式開型求解公式 1. 亞當(dāng)斯顯式求解公式 類似地取類似地取xk4, xk3, xk2, xk1 在在xk4, xk上作上作F(x)=f(x, y(x)的插值多項(xiàng)式，可導(dǎo)出的插值多項(xiàng)式，可導(dǎo)出4 階亞當(dāng)斯顯式求解公式：階亞當(dāng)斯顯式求解公式： (3.6) (3.7) 4階精度階精度 ).,(9),(37),(59),(55 24 443322111 kkkkkkkkkk yxFyxFyxFyxF h yy )()( 720 251 )()( 55 2 )5( hOhyyxyhe kkk 3 線性多步法 v3.1 開型求解公式開型求解公式 2. 亞當(dāng)斯隱式求解公式 取取xk3, xk2

49、, xk1, xk，在，在xk3, xk上作上作F(x) = f(x, y(x) 的插值多項(xiàng)式的插值多項(xiàng)式 用上述方法可導(dǎo)出：用上述方法可導(dǎo)出： (3.8) (3.9) 稱為亞當(dāng)斯隱式求解公式稱為亞當(dāng)斯隱式求解公式. .)( 3 0 3 0 3 i ij j ik ji j F xx xx xL ),(),(5),(19),(9 24 3322111 kkkkkkkkkk yxfyxfyxfyxf h yy )()( 720 19 )()( 55)5( hOhyyxyhe kkk 3 線性多步法 v3.1 開型求解公式開型求解公式 2. 亞當(dāng)斯隱式求解公式 (3.8) (3.9) 稱為亞當(dāng)斯隱

50、式求解公式稱為亞當(dāng)斯隱式求解公式. 注：利用注：利用4階公式階公式(3.6)顯化之：顯化之： (3.10) 稱稱(3.10)為亞當(dāng)斯預(yù)測為亞當(dāng)斯預(yù)測校正系統(tǒng)校正系統(tǒng). ),(),(5),(19),(9 24 3322111 kkkkkkkkkk yxfyxfyxfyxf h yy )()( 720 19 )()( 55)5( hOhyyxyhe kkk ),(),(5),(19),(9 24 9375955 24 3322111 43211 kkkkkkkkkk kkkkkk yxfyxfyxfyxf h yy ffff h yy 3 線性多步法 v3.2 閉型求解系統(tǒng)閉型求解系統(tǒng) 將將xk1

51、, xk擴(kuò)充為擴(kuò)充為xk4, xk，取，取xk4，xk3，xk2，xk1 為節(jié)點(diǎn)，作為節(jié)點(diǎn)，作F(x) = f(x, y(x) 的牛頓前插多項(xiàng)式的牛頓前插多項(xiàng)式. 則則 ).( ! 4 )3)(2)(1( ! 3 )2)(1( ! 2 )1( )()()( )5( 4 3 4 2 44 33 hy tttt F ttt F tt FtF xRxNxF kkkk k k k k k k x x x x k x x kk dxxRdxxNy dxxyxfyy 44 4 )()( )(,( 334 4 3 線性多步法 v3.2 閉型求解系統(tǒng)閉型求解系統(tǒng) 將將xk1, xk擴(kuò)充為擴(kuò)充為xk4, xk，

52、取，取xk4，xk3，xk2，xk1 為節(jié)點(diǎn)，作為節(jié)點(diǎn)，作F(x) = f(x, y(x) 的牛頓前插多項(xiàng)式的牛頓前插多項(xiàng)式. 則則 令令x = xk + (t 4)h 則則 k k k k x x x x kk dxxRdxxNyy 44 )()( 334 4 0 4 3 4 2 443 ! 3 )2)(1( ! 2 )1( )( 4 hdtF ttt F tt FtFdxxN kkkk x x k k .2563 3 4 3 8 3 20 84 4 3 4 2 44 4 3 4 2 44 kkkk kkkk FFFFh FFFFh 3 線性多步法 v3.2 閉型求解系統(tǒng)閉型求解系統(tǒng) 令令x

53、 = xk + (t 4)h 則則 由由 .2563 3 4 )( 4 3 4 2 443 4 kkkk x x FFFFhdxxN k k ,2, 4324 2 434 kkkkkkk FFFFFFF .33 43214 3 kkkkk FFFFF ).,(2),(),(2 3 4 22 3 4 )( 332211 3213 4 kkkkkk kkk x x yxfyxfyxfh FFFhdxxN k k 3 線性多步法 v3.2 閉型求解系統(tǒng)閉型求解系統(tǒng) 令令 (3.11) 稱稱(3.11)為米爾恩求解公式為米爾恩求解公式(Miline). 余項(xiàng)：余項(xiàng)： ).,(2),(),(2 3 4 )( 3322113 4 kkkkkk x x yxfyxfyxfhdxx