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

1、常微分方程數(shù)值解法第1頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.0 基本概念1. 常微分方程的初值問題：稱為具有初值(1.2)的常微分方程. 若f(x,y)在axb, |y|+上連續(xù)，且關(guān)于y滿足Lip條件：常數(shù)L使| f(x, y1) f(x, y2)| L|y1 y2|則初值問題(1.1)(1.2)存在唯一連續(xù)可微解y(x).注：以下總假設(shè)f 滿足Lip條件.第2頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.0 基本概念1. 常微分方程的初值問題：稱為具有初值(1.2)的常微分方程. (1.1)(1.2)等價于微分方程： (1.3)注：一

2、般無初等解(解析解)，即使有形式也復(fù)雜.第3頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.0 基本概念2. 初值問題的數(shù)值解 設(shè)(1.1)(1.2)的解y(x)在節(jié)點(diǎn)xi處的近似解值為 yi y(xi), a x1 x2 xn = b則稱yi (i = 1, 2, , n)為(1.1)(1.2)的數(shù)值解，又稱y(xi)的計(jì)算值.第4頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.0 基本概念3. 數(shù)值方法 兩種轉(zhuǎn)化： 由微分出發(fā)的數(shù)值方法. 由積分 出發(fā)的數(shù)值方法. 計(jì)算方法 步進(jìn)法：從初始條件出發(fā)，逐步求y1, y2, , yn. 又有兩種：單步

3、法，多步法.注：采用等距節(jié)點(diǎn)：第5頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.1 基于數(shù)值微分的求解公式. (1.6)第6頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.1 基于數(shù)值微分的求解公式.1. 前進(jìn)歐拉公式 (1.6)的前半部分為：令 yi+1 = yi + hf(xi, yi) (1.7)其中yi = y(xi) , 則yi+1 y(xi+1)第7頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.1 基于數(shù)值微分的求解公式.1. 前進(jìn)歐拉公式 令 yi+1 = yi + hf(xi, yi) (1.7)其中yi =

4、 y(xi) , 則yi+1 y(xi+1)記 (1.8)則稱(1.7)為前進(jìn)歐拉求解公式. 簡稱為歐拉公式或歐拉法. (1.8)稱為歐拉公式的余項(xiàng)：ei+1(h) = y(xi+1) yi+1 第8頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.1 基于數(shù)值微分的求解公式.2. 后退歐拉公式 (1.6)的后半部分令 yi+1 = yi + hf(xi+1, yi+1) (1.9)其中yi = y(xi), 則yi+1 y(xi+1) 第9頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.1 基于數(shù)值微分的求解公式.2. 后退歐拉公式令 yi+1 =

5、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) (1.10)第10頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.1 基于數(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ū)別：(1.7)為直接計(jì)算公式稱顯式公式

6、.(1.9)為關(guān)于函數(shù)方程稱隱式公式.第11頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.1 基于數(shù)值微分的求解公式.【例1】取h=0.1求解初值問題： (1.11).解： ，xi = ih = 0.1i, (i = 0, 1, 2, , 10) 歐拉法：第12頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.1 基于數(shù)值微分的求解公式.【例1】取h=0.1求解初值問題： (1.11).解： ，xi = ih = 0.1i, (i = 0, 1, 2, , 10) 后退歐拉法： 第13頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言

7、1.1 基于數(shù)值微分的求解公式.注：為避免求解函數(shù)方程，采用顯式與隱式結(jié)合的方法： 此方法稱為 預(yù)測校正系統(tǒng). 求解過程為：第14頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.1 基于數(shù)值微分的求解公式.預(yù)測校正系統(tǒng)：【例2】利用預(yù)測校正系統(tǒng)求解例1.第15頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.1 基于數(shù)值微分的求解公式.預(yù)測校正系統(tǒng)：注：顯式比隱式方便，但有時隱式效果比顯式好.(4介紹).第16頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.2 截?cái)嗾`差定義1.1 稱ek(h) = y(xk) yk為計(jì)算yk的公式

8、第k步的局部截?cái)嗾`差. 注：“局部”是指在計(jì)算第k步時，假定前面yi = y(xi) (i k).而yk y(xk) 歐拉法. 后退歐拉法.一般根據(jù)y(xk)對y(k), y(k)做估計(jì).第17頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.2 截?cái)嗾`差定義1.2 設(shè)ei(h) (i = 1, 2, , k)為求解公式第i步的局部截?cái)嗾`差.稱為該求解公式在點(diǎn)上的整體截?cái)嗾`差.注：局部截?cái)嗾`差ek(h)與yk有關(guān). 整體截?cái)嗾`差Ek(h)與y1, y2, , yk有關(guān).所有ek(h)都與h有關(guān).第18頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.2

9、截?cái)嗾`差定義1.3 若局部截?cái)嗾`差e(h)=O(hp+1)，則稱該求解公式具有p階精度.注：歐拉法具有一階精度.(精度越高越好)第19頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言作業(yè) P208 1，2，3.第20頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.3 基于數(shù)值積分的求解公式 (1.13)若已知y(xk) = yk, 則計(jì)算積分可求出y(xk+1) . 如用矩形公式求積分則有y(xk+1) = y(xk) + hf(xk, yk)令yk+1 = y(xk) + hf(xk, yk)即為歐拉公式. 故歐拉公式又稱矩形法.第21頁，共72頁，2

10、022年，5月20日，20點(diǎn)34分，星期三1 引言1.3 基于數(shù)值積分的求解公式 (1.13)考慮1. 梯形公式記 (1.14)第22頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.3 基于數(shù)值積分的求解公式1. 梯形公式記 (1.14)稱(1.14)為梯形(求解)公式. 簡稱梯形法.第23頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.3 基于數(shù)值積分的求解公式1. 梯形公式梯形(求解)公式, 簡稱梯形法: (1.14)注：梯形公式的余項(xiàng)： 故是二階精度.第24頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1.3 基于數(shù)值積分的求解公式

11、1. 梯形公式 (1.14) 梯形公式為隱式公式.預(yù)測校正系統(tǒng) (1.15)稱(1.15)為改進(jìn)的歐拉公式，也可記為1 引言第25頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.3 基于數(shù)值積分的求解公式1. 梯形公式 (1.14) 可以證明，改進(jìn)歐拉公式也具有二階精度.第26頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.3 基于數(shù)值積分的求解公式【例3】用歐拉法，梯形法以及改進(jìn)歐拉法求解取h=0.1.計(jì)算到x=0.5.解：f(x, y) = xy + 1, a = x0 = 0, b = 0.5, y0 = 1, n = 5(Euler法) 求解

12、公式：yk =yk1+h(xk1yk1+1)= hxk1+(1 h)yk1 + h = 0.1xk1+0.9yk1+0.1 第27頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.3 基于數(shù)值積分的求解公式【例3】用歐拉法，梯形法以及改進(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，得方程第28頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.3 基于數(shù)值積分的求解公式【例3】用歐拉法，梯形法以及

13、改進(jìn)歐拉法求解解：f(x, y) = xy + 1, a = x0 = 0, b = 0.5, y0 = 1, n = 5(改進(jìn)Euler法)求解公式：yk=yk1+h(xk1yk1+1) + xk (yk +h(xkyk+1)+1/2得=0.905yk1+0.045xk1+0.05xk+0.095方程第29頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.3 基于數(shù)值積分的求解公式2. 辛卜生公式 記 (1.17)第30頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.3 基于數(shù)值積分的求解公式2. 辛卜生公式記 (1.17)其余項(xiàng)第31頁，共72頁，

14、2022年，5月20日，20點(diǎn)34分，星期三1 引言1.3 基于數(shù)值積分的求解公式2. 辛卜生公式記 (1.17)將xk1, xk 對分： 調(diào)整下標(biā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)第32頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.3 基于數(shù)值積分的求解公式2. 辛卜生公式記 (1.17) (1.19)稱(1.19)為辛卜生求

15、解公式，其中fi2= f(xi2, y(xi2)，fi1 = f(xi1, y(xi1)，fi = f(xi, y(xi)注： (1.19)的誤差：第33頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三1 引言1.3 基于數(shù)值積分的求解公式2. 辛卜生公式記 (1.17) (1.19)稱(1.19)為辛卜生求解公式，其中fi2= f(xi2, y(xi2)，fi1 = f(xi1, y(xi1)，fi = f(xi, y(xi)注： 隱式(需顯化)多步將在3中討論.第34頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.0 原理 其中K =

16、 f(, y() = y()稱為y在xi1, xi上的平均斜率.歐拉法：改進(jìn)歐拉法：(2.1)第35頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.0 原理 其中K = f(, y() = y()稱為y在xi1, xi上的平均斜率.對(1.17)顯化：辛卜生： (2.4)第36頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.0 原理其中K = f(, y() = y()稱為y在xi1, xi上的平均斜率.設(shè)想：在中多計(jì)算(預(yù)測)幾個點(diǎn)上的值然后可加權(quán)取平均值作為的近似值可能構(gòu)成更高階的公式.一階二階三階

17、第37頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式 (*)其中0 j 1，yi1 +jh是y(xi1 + jh) 的預(yù)測值. 稱(*)為R-K公式注：(2.1)(2.4)分別稱為二階，三階R-K公式. j，j，j為待定系數(shù). 使(*)的階數(shù)盡量高.第38頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式參數(shù)的確定，以m = 2為例. 欲求1，2，2 .原則: 使ei(h) = y(xi) yi的階數(shù)盡可能高第39頁，共72頁，2022

18、年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式展開展開 原則: 使ei(h) = y(xi) yi的階數(shù)盡可能高第40頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式 原則: 使ei(h) = y(xi) yi的階數(shù)盡可能高第41頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式 欲求截?cái)嗾`差ei(h) = y(xi) yi關(guān)于h的階數(shù)盡可能高，應(yīng)使無窮多解，從而有許多2階R-K

19、公式第42頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式應(yīng)使注： 取1= 2= 1/2，2 = 1，即為改進(jìn)歐拉公式.第43頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式應(yīng)使注：取1= 0，2 = 1，2 = 1/2，即為中點(diǎn)公式第44頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式應(yīng)使注：二階R-K公式的截?cái)嗾`差為故為二階方法.相仿可得更高階的R-K公

20、式.第45頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.2 經(jīng)典R-K公式 在4解R-K公式中最重要的是經(jīng)典R-K公式. (2.6)注： (2.6)為4階方法.第46頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.2 經(jīng)典R-K公式 在4解R-K公式中最重要的是經(jīng)典R-K公式. (2.6)注：R-K法對4階以上不一定能提高整數(shù)階.第47頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法2.2 經(jīng)典R-K公式【例4】使用三階，四階R-K法求解初值問題： 的部分計(jì)算

21、值y1，y2，y3，其中h=0.1.解 使用三階R-K法第48頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法【例4】使用三階，四階R-K法求解初值問題： 的部分計(jì)算值y1，y2，y3，其中h=0.1.解 使用三階R-K法第49頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法【例4】使用三階，四階R-K法求解初值問題： 的部分計(jì)算值y1，y2，y3，其中h=0.1.解 使用四階R-K法第50頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法【例4】使用三階，四階R-K法求

22、解初值問題： 的部分計(jì)算值y1，y2，y3，其中h=0.1.解 使用四階R-K法第51頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三2 Runge - Kutta法注 使用R-K法要求具備較好的光滑性，否則效果不如低階的.作業(yè)P209 8 9，10.第52頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三3 線性多步法單步法的優(yōu)點(diǎn)：簡單，計(jì)算yk+1只用yk.缺點(diǎn): 沒有充分利用前面的信息且計(jì)算y(xk+h)較困難回顧Simpson： (1.19)考慮： (3.1)兩種插值求積： 將xk1, xk增加內(nèi)部節(jié)點(diǎn)，改為xk2, xk導(dǎo)出的公式稱為閉型求解公式.線性多步第53頁，

23、共72頁，2022年，5月20日，20點(diǎn)34分，星期三3 線性多步法考慮： (3.1)兩種插值求積： 將xk1, xk增加內(nèi)部節(jié)點(diǎn)，改為xk2, xk導(dǎo)出的公式稱為閉型求解公式.在xk1, xk外增加插值節(jié)點(diǎn)，導(dǎo)出的公式稱為開型求解公式.開型有顯和隱，閉型也有顯和隱.第54頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三3 線性多步法3.1 開型求解公式1. 亞當(dāng)斯顯式求解公式 取節(jié)點(diǎn)xk3, xk2, xk1，在xk3, xk上作F(x) = f(x, y(x) 的插值多項(xiàng)式.第55頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三3 線性多步法3.1 開型求解公式1. 亞

24、當(dāng)斯顯式求解公式 取節(jié)點(diǎn)xk3, xk2, xk1，在xk3, xk上記xki = xk ih, x = xk + th，則 第56頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三3 線性多步法3.1 開型求解公式1. 亞當(dāng)斯顯式求解公式 取節(jié)點(diǎn)xk3, xk2, xk1，記xki = xk ih, x = xk + th，則 代入(3.1)得第57頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三3 線性多步法3.1 開型求解公式1. 亞當(dāng)斯顯式求解公式 取節(jié)點(diǎn)xk3, xk2, xk1，記xki = xk ih, x = xk + th，則 令 (3.4)稱(3.4)為亞

25、當(dāng)斯顯式求解公式(線性多步).第58頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三3 線性多步法3.1 開型求解公式1. 亞當(dāng)斯顯式求解公式 取節(jié)點(diǎn)xk3, xk2, xk1，記xki = xk ih, x = xk + th，則余項(xiàng)： 從而(3.4)具有3階精度. 稱為3階亞當(dāng)斯求解公式.第59頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三3 線性多步法3.1 開型求解公式1. 亞當(dāng)斯顯式求解公式類似地取xk4, xk3, xk2, xk1 在xk4, xk上作F(x)=f(x, y(x)的插值多項(xiàng)式，可導(dǎo)出4階亞當(dāng)斯顯式求解公式： (3.6) (3.7)4階精度第6

26、0頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三3 線性多步法3.1 開型求解公式2. 亞當(dāng)斯隱式求解公式 取xk3, xk2, xk1, xk，在xk3, xk上作F(x) = f(x, y(x) 的插值多項(xiàng)式用上述方法可導(dǎo)出： (3.8) (3.9)稱為亞當(dāng)斯隱式求解公式.第61頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三3 線性多步法3.1 開型求解公式2. 亞當(dāng)斯隱式求解公式 (3.8) (3.9)稱為亞當(dāng)斯隱式求解公式.注：利用4階公式(3.6)顯化之： (3.10) 稱(3.10)為亞當(dāng)斯預(yù)測校正系統(tǒng).第62頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三3 線性多步法3.2 閉型求解系統(tǒng) 將xk1, xk擴(kuò)充為xk4, xk，取xk4，xk3，xk2，xk1為節(jié)點(diǎn)，作F(x) = f(x, y(x) 的牛頓前插多項(xiàng)式. 則 第63頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三3 線性多步法3.2 閉型求解系統(tǒng) 將xk1, xk擴(kuò)充為xk4, xk，取xk4，xk3，xk2，xk1為節(jié)點(diǎn)，作F(x) = f(x, y(x) 的牛頓前插多項(xiàng)式.則 令x = xk + (t 4)h 則 第64頁，共72頁，2022年，5月20日，20點(diǎn)34分，星期三3