常微分方程数值解法

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