常微分方程數(shù)值解法（課堂PPT）

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

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

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

4、分的求解公式基于數(shù)值微分的求解公式. (1.6).(2)()(1)().(2)()(1)(111iiiiiiiiyhxyxyhxyyhxyxyhxy ,)(2)()(2)()()(1111 iiiiiiiiiixxyhxyyhxyxyxyh 1 引言v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式.1. 前進歐拉公式 (1.6)的前半部分為：的前半部分為：令令 yi+1 = yi + hf(xi, yi) (1.7)其中其中yi = y(xi) , 則則yi+1 y(xi+1).(2)(,()()(11iiiiiyhxyxfxyxyh ).(2)(,()()(21iiiiiyhxyxhf

5、xyxy ,)(2)()(2)()()(1111 iiiiiiiiiixxyhxyyhxyxyxyh 1 引言v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式.1. 前進歐拉公式 令令 yi+1 = yi + hf(xi, yi) (1.7)其中其中yi = y(xi) , 則則yi+1 y(xi+1)記記 (1.8)則則稱稱(1.7)為前進歐拉求解公式為前進歐拉求解公式. 簡稱為歐拉公式或歐拉法簡稱為歐拉公式或歐拉法. (1.8)稱為歐拉公式的余項：稱為歐拉公式的余項：ei+1(h) = y(xi+1) yi+1 ).(2)(,()()(21iiiiiyhxyxhfxyxy )(2)(

6、21iiyhhe ).(2)(21iixyhhe 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)()()(111iiiiyhxyxyxyh ).(2)(,()()(2111iiiiiyhxyxhfxyxy ,)(2)()(2)()()(1111 iiiiiiiiiixxyhxyyhxyxyxyh 1 引言v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式.2. 后退歐拉公式令令 yi

7、+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) 余項余項 (1.10).(2)(2),()(,()()(1111111iiiiiiiiiyhyhyxfxyxfhyxyhe 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)為后退歐拉公式為后退歐拉公式(后退

8、歐拉法后退歐拉法). 稱稱(1.10)為后退歐拉法的誤差近似值為后退歐拉法的誤差近似值. 歐拉法與后退歐拉公式的區(qū)別：歐拉法與后退歐拉公式的區(qū)別：(1.7)為直接計算公式稱顯式公式為直接計算公式稱顯式公式.(1.9)為關于函數(shù)方程稱隱式公式為關于函數(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 , 02yxyxyyyxyyxf2),( . 9,.,2 , 1 , 0)2(1 iy

9、xyhyyiiiii1 引言v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式.【例例1】取取h=0.1求解初值問題：求解初值問題： (1.11).解：解： ，xi = ih = 0.1 i, (i = 0, 1, 2, , 10) 后退歐拉法：后退歐拉法： 1)0(1 , 02yxyxyyyxyyxf2),( )2(1111 iiiiiyxyhyy. 9,.,2 , 1 , 0)1(2)1(8121 ihhxhyyyiiii1 引言v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式.注：為避免求解函數(shù)方程，采用顯式與隱式結合的方注：為避免求解函數(shù)方程，采用顯式與隱式結合的方法：法：

10、此方法稱為此方法稱為 預測預測校正系統(tǒng)校正系統(tǒng). 求解過程為：求解過程為：校校正正值值隱隱式式預預測測值值顯顯式式)()(),(),(1111 iiiiiiiiyxhfyyyxhfyy.)()()(22110nnyyyyyyy1 引言v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式.預測預測校正系統(tǒng)：校正系統(tǒng)：【例例2】利用預測利用預測校正系統(tǒng)求解例校正系統(tǒng)求解例1.校校正正值值隱隱式式預預測測值值顯顯式式)()(),(),(1111 iiiiiiiiyxhfyyyxhfyy 1)0(1 , 02yxyxyy. 9,.,2 , 1 , 0)2()2(11111 iyxyhyyyxyhyy

11、iiiiiiiiii1 引言v1.1 基于數(shù)值微分的求解公式基于數(shù)值微分的求解公式.預測預測校正系統(tǒng)：校正系統(tǒng)：注：顯式比隱式方便，但有時隱式效果比顯式好注：顯式比隱式方便，但有時隱式效果比顯式好.(4介紹介紹).校校正正值值隱隱式式預預測測值值顯顯式式)()(),(),(1111 iiiiiiiiyxhfyyyxhfyy1 引言v1.2 截斷誤差截斷誤差定義定義1.1 稱稱ek(h) = y(xk) yk為計算為計算yk的公式第的公式第k步的局步的局部截斷誤差部截斷誤差. 注：注：“局部局部”是指在計算第是指在計算第k步時，假定前面步時，假定前面yi = y(xi) (i k).而而yk y

12、(xk) 歐拉法歐拉法. 后退歐拉法后退歐拉法.一般根據(jù)一般根據(jù)y(xk)對對y( k), y( k)做估計做估計.)(2)(21kkyhhe ).(2)(21iiyhhe 1 引言v1.2 截斷誤差截斷誤差定義定義1.2 設設ei(h) (i = 1, 2, , k)為求解公式第為求解公式第i步的局部步的局部截斷誤差截斷誤差.稱稱為該求解公式在點上的整體截斷誤差為該求解公式在點上的整體截斷誤差.注：局部截斷誤差注：局部截斷誤差ek(h)與與yk有關有關. 整體截斷誤差整體截斷誤差Ek(h)與與y1, y2, , yk有關有關.所有所有ek(h)都與都與h有關有關. kiikhehE1)()(

13、1 引言v1.2 截斷誤差截斷誤差定義定義1.3 若局部截斷誤差若局部截斷誤差e(h)=O(hp+1)，則稱該求解公，則稱該求解公式具有式具有p階精度階精度.注：歐拉法具有一階精度注：歐拉法具有一階精度.(精度越高越好精度越高越好)1 引言作業(yè)作業(yè) P208 1，2，3.1 引言v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 (1.13)若已知若已知y(xk) = yk, 則計算積分可求出則計算積分可求出y(xk+1) . 如用矩形公式求積分如用矩形公式求積分則有則有y(xk+1) = y(xk) + hf(xk, yk)令令yk+1 = y(xk) + hf(xk, yk)即為歐拉公

14、式即為歐拉公式. 故歐拉公式又故歐拉公式又稱矩形法稱矩形法.).()()()(,(111kkxxxxxyxydxxydxxyxfkkkk .)(,()()(11 kkxxkkdxxyxfxyxy).,()(,(1kkxxxxhfdxxyxfkk 1 引言v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式 (1.13)考慮考慮1. 梯形公式記記 (1.14).)(,()()(11 kkxxkkdxxyxfxyxy 11)()(,(kkkkxxxxdxxFdxxyxf).(,()(,(2)()(,()(,(2)(,(111111 kkkkkkkkkkxxxyxfxyxfhyxyxyxfxyxf

15、hdxxyxfkk),(),(2111 kkkkkkyxfyxfhyy1 引言v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式1. 梯形公式記記 (1.14)稱稱(1.14)為梯形為梯形(求解求解)公式公式. 簡稱梯形法簡稱梯形法.).(,()(,(2)()(,()(,(2)(,(111111 kkkkkkkkkkxxxyxfxyxfhyxyxyxfxyxfhdxxyxfkk),(),(2111 kkkkkkyxfyxfhyy1 引言v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式1. 梯形公式梯形梯形(求解求解)公式公式, 簡稱梯形法簡稱梯形法: (1.14)注：梯形公式的余項：

16、注：梯形公式的余項： 故是二階精度故是二階精度.),(),(2111 kkkkkkyxfyxfhyy,)(12)()(1331 iiiiixxhOhyhe 3)(12)(abfRT v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式1. 梯形公式 (1.14) 梯形公式為隱式公式梯形公式為隱式公式.預測預測校正系統(tǒng)校正系統(tǒng) (1.15)稱稱(1.15)為改進的歐拉公式，也可記為為改進的歐拉公式，也可記為1 引言),(),(2111 kkkkkkyxfyxfhyy. 1,.,2 , 1 , 0)()(),(),(1111 niyxhfyyyxhfyyiiiiiiii梯形梯形校正校正歐拉歐拉預

17、測預測),(,(),(211iiiiiiiiyxhfyxfyxfhyy 1 引言v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式1. 梯形公式 (1.14) 可以證明，改進歐拉公式也具有二階精度可以證明，改進歐拉公式也具有二階精度.),(),(2111 kkkkkkyxfyxfhyy1 引言v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式【例例3】用歐拉法，梯形法以及改進歐拉法求解用歐拉法，梯形法以及改進歐拉法求解取取h=0.1.計算到計算到x=0.5.解：解：f(x, y) = xy + 1, a = x0 = 0, b = 0.5, y0 = 1, n = 5(Euler法法)

18、 求解公式：求解公式：yk =yk1+h(xk1yk1+1)= hxk1+(1 h)yk1 + h = 0.1xk1+0.9yk1+0.1 1)0(1yxyy1 引言v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式【例例3】用歐拉法，梯形法以及改進歐拉法求解用歐拉法，梯形法以及改進歐拉法求解解：解：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(1yxyy1 . 22 . 0)( 1 . 09 . 122)()

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

20、 引言v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式2. 辛卜生公式 記記 (1.17)(,()2(,2(4)(,(6)(,(11111kkkkkkkkxxxyxfxxyxxfxyxfhdxxyxfkk )(,()2(,2(4)(,(61111hxyhxfhxyhxfxyxfhkkkkkk )(,()2(,2(4)(,(611111hxyhxfhxyhxfxyxfhyykkkkkkkk 1 引言v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式2. 辛卜生公式記記 (1.17)其余項其余項)(,()2(,2(4)(,(611111hxyhxfhxyhxfxyxfhyykkkkkkk

21、k )(,(16180)()4(5kkkyfhhe ,)()(1618015)5(5iiiixxhOyh 1 引言v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式2. 辛卜生公式記記 (1.17)將將xk1, xk 對分：對分： 調(diào)整下標為調(diào)整下標為xi2, xi ：xi2 = xk1, xi1 = xk1+h1, xi = xk1+2h1= xk則則(1.17)化為化為 (1.19)稱稱(1.19)為辛卜生求解公式，其中為辛卜生求解公式，其中fk2= f(xk2, y(xk2)，fk1 = f(xk1, y(xk1)，fk = f(xk, y(xk)(,()2(,2(4)(,(6111

22、11hxyhxfhxyhxfxyxfhyykkkkkkkk .21hh 431212iiiiifffhyy 1 引言v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式2. 辛卜生公式記記 (1.17) (1.19)稱稱(1.19)為辛卜生求解公式，其中為辛卜生求解公式，其中fi2= f(xi2, y(xi2)，fi1 = f(xi1, y(xi1)，fi = f(xi, y(xi)注：注： (1.19)的誤差：的誤差：)(,()2(,2(4)(,(611111hxyhxfhxyhxfxyxfhyykkkkkkkk 431212iiiiifffhyy );()(90)(16180)2()(5

23、1)(51)(511hOyhyhhekkkkk 1 引言v1.3 基于數(shù)值積分的求解公式基于數(shù)值積分的求解公式2. 辛卜生公式記記 (1.17) (1.19)稱稱(1.19)為辛卜生求解公式，其中為辛卜生求解公式，其中fi2= f(xi2, y(xi2)，fi1 = f(xi1, y(xi1)，fi = f(xi, y(xi)注：注： 隱式隱式(需顯化需顯化)多步多步將在將在3中討論中討論.)(,()2(,2(4)(,(611111hxyhxfhxyhxfxyxfhyykkkkkkkk 431212iiiiifffhyy 2 Runge - Kutta法v2.0 原理原理 其中其中K = f(

24、 , y( ) = y( )稱為稱為y在在xi1, xi上的平均斜率上的平均斜率.歐拉法：歐拉法：改進歐拉法：改進歐拉法：(2.1)(,()()(,()()(111 yhfxydxxyxfxyxyixxiiii hKxyi )(1 ),(11111iiiiyxfKhKyy ),(),()2(112111211hKyxfKyxfKKKhyyiiiiii2 Runge - Kutta法v2.0 原理原理 其中其中K = f( , y( ) = y( )稱為稱為y在在xi1, xi上的平均斜率上的平均斜率.對對(1.17)顯化：顯化：辛卜生：辛卜生： (2.4)(,()()(,()()(111 yh

25、fxydxxyxfxyxyixxiiii hKxyi )(1 )2(,()2,2(),(641211311121113211KKhyhxfKKhyhxfKyxfKKKKhyyiiiiiiii2 Runge - Kutta法v2.0 原理原理其中其中K = f( , y( ) = y( )稱為稱為y在在xi1, xi上的平均斜率上的平均斜率.設想：在中多計算設想：在中多計算(預測預測)幾個點上的值然后可加權取幾個點上的值然后可加權取平均值作為的近似值可能構成更高階的公式平均值作為的近似值可能構成更高階的公式. hKxyxyii )()(1 ),(11111iiiiyxfKhKyy ),(),()

26、2(112111211hKyxfKyxfKKKhyyiiiiii )2(,()2,2(),(641211311121113211KKhyhxfKKhyhxfKyxfKKKKhyyiiiiiiii一階一階二階二階三階三階2 Runge - Kutta法v2.1 Runge - Kutta公式公式 (*)其中其中0 j 1，yi1 + jh是是y(xi1 + jh) 的預測值的預測值. 稱稱(*)為為R-K公式公式注：注：(2.1)(2.4)分別稱為二階，三階分別稱為二階，三階R-K公式公式. j， j， j為待定系數(shù)為待定系數(shù). 使使(*)的階數(shù)盡量高的階數(shù)盡量高. ),(),(),()(112

27、121211122111hyhxfKhyhxfKyxfKKKKhyymimimiiiimmii 2 Runge - Kutta法v2.1 Runge - Kutta公式公式參數(shù)的確定，以參數(shù)的確定，以m = 2為例為例. 欲求欲求 1， 2， 2 . ),()(,(),(121212121211122111hKyhxfhxyhxfKyxfKKKhyyiiiiiiii 原則原則: 使使ei(h) = y(xi) yi的階數(shù)盡可能高的階數(shù)盡可能高2 Runge - Kutta法v2.1 Runge - Kutta公式公式展開展開展開展開 ),()(,(),(121212121211122111hK

28、yhxfhxyhxfKyxfKKKhyyiiiiiiii 原則原則: 使使ei(h) = y(xi) yi的階數(shù)盡可能高的階數(shù)盡可能高).()(2)()()()(312111hOxyhxyhxyhxyxyiiiii )(),(),(211122112hOyxfyhKxhyxfKiiii )(),(),(),(211111211hOyxfyKyxfxhyxfiiiiii 2 Runge - Kutta法v2.1 Runge - Kutta公式公式 ),()(,(),(121212121211122111hKyhxfhxyhxfKyxfKKKhyyiiiiiiii 原則原則: 使使ei(h) =

29、y(xi) yi的階數(shù)盡可能高的階數(shù)盡可能高).()(2)()()()(312111hOxyhxyhxyhxyxyiiiii )(),(),(),(2111112112hOyxfyKyxfxhyxfKiiiiii . )()()()()(),(),(),(),(3122212113111112221121111hOxyhxyhxyhOyxfyKyxfxhyxhfyxhfyyiiiiiiiiiiiii 2 Runge - Kutta法v2.1 Runge - Kutta公式公式 欲求截斷誤差欲求截斷誤差ei(h) = y(xi) yi關于關于h的階數(shù)盡可能高，的階數(shù)盡可能高，應使應使 ),()(

30、,(),(121212121211122111hKyhxfhxyhxfKyxfKKKhyyiiiiiiii 無窮多解，從而有許多無窮多解，從而有許多2階階R-K公式公式).()(2)()()()(312111hOxyhxyhxyhxyxyiiiii ).()()()(312221211hOxyhxyhxyyiiii 2112221 2 Runge - Kutta法v2.1 Runge - Kutta公式公式應使應使注：注： 取取 1= 2= 1/2， 2 = 1，即為改進歐拉公式，即為改進歐拉公式. ),()(,(),(121212121211122111hKyhxfhxyhxfKyxfKKK

31、hyyiiiiiiii 2112221 ),(),()2(112111211hKyxfKyxfKKKhyyiiiiii2 Runge - Kutta法v2.1 Runge - Kutta公式公式應使應使注：注：取取 1= 0， 2 = 1， 2 = 1/2，即為中點公式，即為中點公式 ),()(,(),(121212121211122111hKyhxfhxyhxfKyxfKKKhyyiiiiiiii 2112221 )2,2(),(111211121KhyhxfKyxfKhKyyiiiiii2 Runge - Kutta法v2.1 Runge - Kutta公式公式應使應使注：二階注：二階R-

32、K公式的截斷誤差為故為二階方法公式的截斷誤差為故為二階方法.相仿相仿可得更高階的可得更高階的R-K公式公式. ),()(,(),(121212121211122111hKyhxfhxyhxfKyxfKKKhyyiiiiiiii 2112221 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(),(22631142113111211143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyyiiiiiiiiii2

33、 Runge - Kutta法v2.2 經(jīng)典經(jīng)典R-K公式公式 在在4解解R-K公式中最重要的是經(jīng)典公式中最重要的是經(jīng)典R-K公式公式. (2.6)注：注：R-K法對法對4階以上不一定能提高整數(shù)階階以上不一定能提高整數(shù)階. ),()2,2()2,2(),(22631142113111211143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyyiiiiiiiiii2 Runge - Kutta法v2.2 經(jīng)典經(jīng)典R-K公式公式【例例4】使用三階，四階使用三階，四階R-K法求解初值問題：法求解初值問題： 的部分計算值的部分計算值y1，y2，y3，其中，其中h=0.1.解解 使

34、用三階使用三階R-K法法 1)0(5 . 0 , 02yxyy 2121121132111112211113211)2()2(,()2()2,2(),(64KKhyKKhyhxfKKhyKhyhxfKyyxfKKKKhyyiiiiiiiiiii2 Runge - Kutta法【例例4】使用三階，四階使用三階，四階R-K法求解初值問題：法求解初值問題： 的部分計算值的部分計算值y1，y2，y3，其中，其中h=0.1.解解 使用三階使用三階R-K法法 1)0(5 . 0 , 02yxyy 2121321122113211)2()2(64KKhyKKhyKyKKKKhyyiiiii2 Runge -

35、 Kutta法【例例4】使用三階，四階使用三階，四階R-K法求解初值問題：法求解初值問題： 的部分計算值的部分計算值y1，y2，y3，其中，其中h=0.1.解解 使用四階使用四階R-K法法 1)0(5 . 0 , 02yxyy 2313114221211321111122111143211)(),()2()2,2()2()2,2()(),(226hKyhKyhxfKKhyKhyhxfKKhyKhyhxfKyyxfKKKKKhyyiiiiiiiiiiiiii2 Runge - Kutta法【例例4】使用三階，四階使用三階，四階R-K法求解初值問題：法求解初值問題： 的部分計算值的部分計算值y1，

36、y2，y3，其中，其中h=0.1.解解 使用四階使用四階R-K法法 1)0(5 . 0 , 02yxyy 23142213211221143211)()2()2()(226hKyKKhyKKhyKyKKKKKhyyiiiiii2 Runge - Kutta法注注 使用使用R-K法要求具備較好的光滑性，否則效果不如法要求具備較好的光滑性，否則效果不如低階的低階的.作業(yè)作業(yè)P209 8 9，10.3 線性多步法單步法的優(yōu)點：簡單，計算單步法的優(yōu)點：簡單，計算yk+1只用只用yk.缺點缺點: 沒有充分利用前面的信息且計算沒有充分利用前面的信息且計算y(xk+ h)較困難較困難回顧回顧Simpson：

37、 (1.19)考慮：考慮： (3.1)兩種插值求積：兩種插值求積： 將將xk1, xk增加內(nèi)部節(jié)點，改為增加內(nèi)部節(jié)點，改為xk2, xk導出的公式導出的公式稱為閉型求解公式稱為閉型求解公式.431212kkkkkfffhyy 線性多步線性多步.)(,()(11 kkxxkkdxxyxfyxy3 線性多步法考慮：考慮： (3.1)兩種插值求積：兩種插值求積： 將將xk1, xk增加內(nèi)部節(jié)點，改為增加內(nèi)部節(jié)點，改為xk2, xk導出的公式導出的公式稱為閉型求解公式稱為閉型求解公式.在在xk1, xk外增加插值節(jié)點，導出的公式稱為開型外增加插值節(jié)點，導出的公式稱為開型求解公式求解公式.開型有顯和隱，

38、閉型也有顯和隱開型有顯和隱，閉型也有顯和隱.)(,()(11 kkxxkkdxxyxfyxy3 線性多步法v3.1 開型求解公式開型求解公式1. 亞當斯顯式求解公式 取節(jié)點取節(jié)點xk3, xk2, xk1，在，在xk3, xk上作上作F(x) = f(x, y(x) 的插值多項式的插值多項式. 3132313212321231131213222)(! 3)()()()()()()()()()()()()(iikkkkkkkkkkkkkkkkkkkkkkxxFxFxxxxxxxxxFxxxxxxxxxFxxxxxxxxxRxLxF 3 線性多步法v3.1 開型求解公式開型求解公式1. 亞當斯顯式

39、求解公式 取節(jié)點取節(jié)點xk3, xk2, xk1，在，在xk3, xk上上記記xki = xk ih, x = xk + th，則，則 3132313212321231131213222)(! 3)()()()()()()()()()()()()(iikkkkkkkkkkkkkkkkkkkkkkxxFxFxxxxxxxxxFxxxxxxxxxFxxxxxxxxxRxLxF .2)2)(1()3)(1(2)3)(2()()(32122 kkkFttFttFtttLxL3 線性多步法v3.1 開型求解公式開型求解公式1. 亞當斯顯式求解公式 取節(jié)點取節(jié)點xk3, xk2, xk1，記，記xki =

40、 xk ih, x = xk + th，則，則 代入代入(3.1)得得.2)2)(1()3)(1(2)3)(2()()(32122 kkkFttFttFtttLxL.51623122)2)(1()3)(1(2)3)(2()()()(321130120110110121211 kkkkkkkkkxxkkFFFhyFdtttFdtttFdttthydttLhydxxLyxykk3 線性多步法v3.1 開型求解公式開型求解公式1. 亞當斯顯式求解公式 取節(jié)點取節(jié)點xk3, xk2, xk1，記，記xki = xk ih, x = xk + th，則，則 令令 (3.4)稱稱(3.4)為亞當斯顯式求解

41、公式為亞當斯顯式求解公式(線性多步線性多步).5162312)(3211 kkkkkFFFhyxy).,(5),(16),(23123322111 kkkkkkkkyxFyxFyxFhyy3 線性多步法v3.1 開型求解公式開型求解公式1. 亞當斯顯式求解公式 取節(jié)點取節(jié)點xk3, xk2, xk1，記，記xki = xk ih, x = xk + th，則，則余項：余項： 從而從而(3.4)具有具有3階精度階精度. 稱為稱為3階亞當斯求解公式階亞當斯求解公式.)1(! 3)()(! 3)()(331)4(313htyxxFxRiiik .)(83)3)(2)(1()(6)1()(6)(4)4

42、(011)4(40131)4(431hydttttyhdttyhdxxRixxkk 3 線性多步法v3.1 開型求解公式開型求解公式1. 亞當斯顯式求解公式類似地取類似地取xk4, xk3, xk2, xk1 在在xk4, xk上作上作F(x)=f(x, y(x)的插值多項式，可導出的插值多項式，可導出4階亞當斯顯式求解公式：階亞當斯顯式求解公式： (3.6) (3.7)4階精度階精度).,(9),(37),(59),(5524443322111 kkkkkkkkkkyxFyxFyxFyxFhyy)()(720251)()(552)5(hOhyyxyhekkk 3 線性多步法v3.1 開型求解

43、公式開型求解公式2. 亞當斯隱式求解公式 取取xk3, xk2, xk1, xk，在，在xk3, xk上作上作F(x) = f(x, y(x) 的插值多項式的插值多項式用上述方法可導出：用上述方法可導出： (3.8) (3.9)稱為亞當斯隱式求解公式稱為亞當斯隱式求解公式.)(30303 iijjikjijFxxxxxL),(),(5),(19),(9243322111 kkkkkkkkkkyxfyxfyxfyxfhyy)()(72019)()(55)5(hOhyyxyhekkk 3 線性多步法v3.1 開型求解公式開型求解公式2. 亞當斯隱式求解公式 (3.8) (3.9)稱為亞當斯隱式求解

44、公式稱為亞當斯隱式求解公式.注：利用注：利用4階公式階公式(3.6)顯化之：顯化之： (3.10) 稱稱(3.10)為亞當斯預測為亞當斯預測校正系統(tǒng)校正系統(tǒng).),(),(5),(19),(9243322111 kkkkkkkkkkyxfyxfyxfyxfhyy)()(72019)()(55)5(hOhyyxyhekkk ),(),(5),(19),(924937595524332211143211 kkkkkkkkkkkkkkkkyxfyxfyxfyxfhyyffffhyy3 線性多步法v3.2 閉型求解系統(tǒng)閉型求解系統(tǒng) 將將xk1, xk擴充為擴充為xk4, xk，取，取xk4，xk3，xk

45、2，xk1為節(jié)點，作為節(jié)點，作F(x) = f(x, y(x) 的牛頓前插多項式的牛頓前插多項式. 則則 ).(! 4)3)(2)(1(! 3)2)(1(! 2)1()()()()5(43424433 hyttttFtttFttFtFxRxNxFkkkk kkkkkkxxxxkxxkkdxxRdxxNydxxyxfyy444)()()(,(33443 線性多步法v3.2 閉型求解系統(tǒng)閉型求解系統(tǒng) 將將xk1, xk擴充為擴充為xk4, xk，取，取xk4，xk3，xk2，xk1為節(jié)點，作為節(jié)點，作F(x) = f(x, y(x) 的牛頓前插多項式的牛頓前插多項式.則則 令令x = xk + (

46、t 4)h 則則 kkkkxxxxkkdxxRdxxNyy44)()(334 404342443! 3)2)(1(! 2)1()(4hdtFtttFttFtFdxxNkkkkxxkk.2563343832084434244434244 kkkkkkkkFFFFhFFFFh3 線性多步法v3.2 閉型求解系統(tǒng)閉型求解系統(tǒng)令令x = xk + (t 4)h 則則由由.256334)(43424434 kkkkxxFFFFhdxxNkk,2,43242434 kkkkkkkFFFFFFF.33432143 kkkkkFFFFF).,(2),(),(2342234)(33221132134 kkkkkkkkkxxyxfyxfyxfhFFFhdxxNkk3 線性多步法v3.2 閉型求解系統(tǒng)閉型求解系統(tǒng)令令 (3.11)稱稱(3.11)為米爾恩求解公式為米爾恩求解公式(Miline).余項：余項：).,(2),(),(234)(33221134 kkkkkkxxyxfyxfyxfhdxxNkk).,(2),(),(2343322114 kkkkkkkkyxfyxfyxfhyy).3(2)2()(2(34)