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

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

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

4、 y(xi) , 則yi+1 y(xi+1)記 (1.8)則稱(1.7)為前進歐拉求解公式. 簡稱為歐拉公式或歐拉法. (1.8)稱為歐拉公式的余項：ei+1(h) = y(xi+1) yi+1 第8頁，共72頁，2022年，5月20日，20點34分，星期三1 引言1.1 基于數值微分的求解公式.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點34分，星期三1 引言1.1 基于數值微分的求解公式.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) 余項 (1.10)第10頁，共72頁，2022年，5月20日，20點34分，星期三1 引言1.1 基于數值微分的求解公式.2. 后退歐拉公式令 yi+1 = yi + hf(xi+1, yi+1) (1.9)其中yi = y(xi), 則yi+1 y(xi+1) 注： 稱(1.9)為后退歐拉公式(后退歐拉法). 稱(1.10)為后退歐拉法的誤差近似值. 歐拉法與后退歐拉公式的區別：(1.7)為直接計算公式稱顯式公式

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

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

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

9、截斷誤差定義1.3 若局部截斷誤差e(h)=O(hp+1)，則稱該求解公式具有p階精度.注：歐拉法具有一階精度.(精度越高越好)第19頁，共72頁，2022年，5月20日，20點34分，星期三1 引言作業 P208 1，2，3.第20頁，共72頁，2022年，5月20日，20點34分，星期三1 引言1.3 基于數值積分的求解公式 (1.13)若已知y(xk) = yk, 則計算積分可求出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點34分，星期三1 引言1.3 基于數值積分的求解公式 (1.13)考慮1. 梯形公式記 (1.14)第22頁，共72頁，2022年，5月20日，20點34分，星期三1 引言1.3 基于數值積分的求解公式1. 梯形公式記 (1.14)稱(1.14)為梯形(求解)公式. 簡稱梯形法.第23頁，共72頁，2022年，5月20日，20點34分，星期三1 引言1.3 基于數值積分的求解公式1. 梯形公式梯形(求解)公式, 簡稱梯形法: (1.14)注：梯形公式的余項： 故是二階精度.第24頁，共72頁，2022年，5月20日，20點34分，星期三1.3 基于數值積分的求解公式

11、1. 梯形公式 (1.14) 梯形公式為隱式公式.預測校正系統 (1.15)稱(1.15)為改進的歐拉公式，也可記為1 引言第25頁，共72頁，2022年，5月20日，20點34分，星期三1 引言1.3 基于數值積分的求解公式1. 梯形公式 (1.14) 可以證明，改進歐拉公式也具有二階精度.第26頁，共72頁，2022年，5月20日，20點34分，星期三1 引言1.3 基于數值積分的求解公式【例3】用歐拉法，梯形法以及改進歐拉法求解取h=0.1.計算到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點34分，星期三1 引言1.3 基于數值積分的求解公式【例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，得方程第28頁，共72頁，2022年，5月20日，20點34分，星期三1 引言1.3 基于數值積分的求解公式【例3】用歐拉法，梯形法以及

13、改進歐拉法求解解：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方程第29頁，共72頁，2022年，5月20日，20點34分，星期三1 引言1.3 基于數值積分的求解公式2. 辛卜生公式 記 (1.17)第30頁，共72頁，2022年，5月20日，20點34分，星期三1 引言1.3 基于數值積分的求解公式2. 辛卜生公式記 (1.17)其余項第31頁，共72頁，

14、2022年，5月20日，20點34分，星期三1 引言1.3 基于數值積分的求解公式2. 辛卜生公式記 (1.17)將xk1, xk 對分： 調整下標為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點34分，星期三1 引言1.3 基于數值積分的求解公式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點34分，星期三1 引言1.3 基于數值積分的求解公式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點34分，星期三2 Runge - Kutta法2.0 原理 其中K =

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

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

18、年，5月20日，20點34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式展開展開 原則: 使ei(h) = y(xi) yi的階數盡可能高第40頁，共72頁，2022年，5月20日，20點34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式 原則: 使ei(h) = y(xi) yi的階數盡可能高第41頁，共72頁，2022年，5月20日，20點34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式 欲求截斷誤差ei(h) = y(xi) yi關于h的階數盡可能高，應使無窮多解，從而有許多2階R-K

19、公式第42頁，共72頁，2022年，5月20日，20點34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式應使注： 取1= 2= 1/2，2 = 1，即為改進歐拉公式.第43頁，共72頁，2022年，5月20日，20點34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式應使注：取1= 0，2 = 1，2 = 1/2，即為中點公式第44頁，共72頁，2022年，5月20日，20點34分，星期三2 Runge - Kutta法2.1 Runge - Kutta公式應使注：二階R-K公式的截斷誤差為故為二階方法.相仿可得更高階的R-K公

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

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

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

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

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

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

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