1 第八章常微分方程初值问题的数值解法

2 第八章 常微分方程数值解法 8.1 引言 ( 基本求解公式 )8.1 引言 ( 基本求解公式 ) 8.2 Runge-Kutta 法8.2 Runge-Kutta 法 8.3 微分方程组和高阶方程解法简介8.3 微分方程组和高阶方程解法简介

3 本章要点：本章主要研究基于微积分数值解法的常 微分方程数值解，主要方法有 线性单步法中的 Euler 方法、 Simpson 方法、 Runge-Kutta 方法 高阶微分方程和微分方程组的数值解法

4 本章应用题 : 驱逐舰在浓雾中搜索潜艇，其时发现潜艇在 3 英里 的海面上，但潜艇立即下潜，驱逐舰速度两倍于潜 艇，且已知潜艇下潜后即以全速朝某一未知方向直 线前进，问驱逐舰应采取什么路线才能保证它会开 过潜艇的上方以投放深水炸弹？ 提示 取极坐标，并以发现潜艇时潜艇的位置为原点 ——— 反潜

5 8.1 引言 ( 基本求解公式 ) 在工程和科学技术的实际问题中，常需要求解微分方程 只有简单的和典型的微分方程可以求出解析解 而在实际问题中的微分方程往往无法求出解析解 在高等数学中我们见过以下常微分方程： (1) (2)

(3) (1),(2) 式称为初值问题,(3) 式称为边值问题 (4) 另外, 在实际应用中还经常需要求解常微分方程组 : 本课程主要研究问题 (1) 的数值解法, 对 (2)~(4) 只作简单介绍 我们首先介绍初值问题 (1) 的解存在的条件

7 定理 1. 对于问题 (1), 要求它的数值解

(1) 从 (1) 的表达式 可以看出, 求它的数值解的关键在于 而数值微分或数值积分问题我们都已经学习过

9 一、基于数值微分的常微分方程数值解法 (1) 对于初值问题 (1) 在下列子区间上分别应用两点数值微分公式 为了讨论方便, 假设以下节点为等距节点

(5) ( 一 ) Euler 公式

11 由 (5) 式每组的前一半可得 (6) (7) 记 其中 (6) 和 (7) 式称为求解初值问题 (1) 的 ( 前进 )Euler 公式和误差项

12 由 (5) 式每组的后一半可得 记 其中 (8) (9) (8) 和 (9) 式称为求解初值问题 (1) 的后退 Euler 公式和误差项

13 从 (6) 或 (8) 式不难看出, 这种类型的方法称为单步格式或单步法 Euler 方法的几何体现 : 前进 Euler 公式 后退 Euler 公式

14 Euler1.m 例 1. 解:解: 由前进 Euler 公式

15 得 依此类推, 有

16 由于后退 Euler 公式是隐形公式, 计算例 1 将很麻烦 事实上大多数情况下用后退 Euler 公式都较困难 就可得到新的 Euler 公式 (10) 此方法称为预测 — 校正系统

17 用 Euler 公式的预测 —— 校正系统求解例 1.例 1 例 2. 解:解: 由 (10) 式, 有 Euler1.m

18 依此类推, 得 比较不同的结果

19 ( 二 ) 常微分方程数值解的截断误差 评价一个微分方程求解公式的标准当然是其精度 而在求解公式 中 误差项

20 定义 1. 因为一般情况下, 求解公式的每一步都存在误差, 因此有 定义 2. 定义 3.

21

22 Euler 公式的局部截断误差为 具有 1 阶精度 后退 Euler 公式的局部截断误差为 也具有 1 阶精度 显然一个求解公式的精度越高, 计算解的精确性也就越好 从前面的分析可知,Euler 法的精度并不算高 因此有必要找寻精度更高的求解公式

23 二、基于数值积分的常微分方程数值解法 (1) 对于初值问题 (11)

24 矩形求积公式 梯形求积公式, 误差为 Simpson 求积公式, 误差为 将以上求积公式代入 (11) 式, 并加以处理就 可得到相对应的求解公式

25 ( 一 ) 矩形求解公式 由 可得 令 (12) (12) 式称为矩形公式 ( 矩形法 ) 实际上就是 Euler 求解公式

26 ( 二 ) 梯形求解公式 由 可得 令 (13) 称 (13) 式为梯形求解公式 ( 梯形法 ) 注意 :(13) 式是隐形公式

27 则梯形公式第 k 步的截断误差为 显然梯形法具有二阶精度 由于梯形公式为隐形公式, 一般情况下不易显化

(14) 以上公式称为改进的 Euler 求解公式 ( 改进 Euler 法 ), 即 (15)

29 例 3. 用 Euler 公式、梯形公式和改进 Euler 公式求 解初值问题, 并比较结果的精度 解:解: (1)Euler 公式

30 (2) 梯形公式

31

32 (3) 改进 Euler 公式 改进 Euler 公式 x y 使用 MATLAB 软件 Euler2.m 结果为

Euler 公式 梯形公式 改进 Euler 公式 结果比较 Euler 法的精度不如梯形公式

34

35 ( 三 ) Simpson 求解公式 将 Simpson 求积公式 代入（ 11 ） 简化后, 得 (16)

36 由 Simpson 求积公式的误差 可以近似得到 (16) 式的截断误差为 仔细分析 (16) 式 如何求 ? 求积公式 (16) 的精度为 4 阶

37 考虑将 (16) 式改为下面的形式 (17) (18) (17) 式称为 Simpson 求解公式,(18) 为相应的截断误差项 (17) 式是一个隐式求解公式 这种形式称为多步法, 在本课程中将不予介绍

Runge-Kutta 法 考虑改进 Euler 法 如果将其改成 (1)

39 改进 Euler 法是由梯形公式和 Euler 公式复合而成 梯形公式具有 2 阶精度 形如 (1) 式的求解公式称为二阶 Runge-Kutta 法 同样可以证明, 改进 Euler 法也具有 2 阶精度 对于 Simpson 求解公式： 这是隐式多步法 选取适当的显化方法, 可得类似 (1) 的高阶 Runge-Kutta 方法 以下使用中值定理进行推导

40 一、 Runge-Kutta 方法的导出 对于常微分方程的边值问题 的解 即 (3) 为了同学们课后复习的方便, 以下的内容将 k 写成 n.

(3) 引入记号 就可得到相应的 Runge-Kutta 方法 (4) 即 (4) 式

42 二、低阶 Runge-Kutta 方法 如下图 即 则 (4) 式化为 即 Euler 方法 Euler 方法也称为一阶 Runge-Kutta 方法 由于 ----(5)

43 ( 由 (5) 式 ) 令 则 (4) 式化为

(6) 和 (1) 式一致, 即改进 Euler 公式, 也称为二阶 Runge-Kutta 法 (1) 式

45 三、高阶 Runge-Kutta 方法 未知

46 令 令 参照 Simpson 求解公式

47 取 则 (4) 式化为 (4) 式 (7) (7) 式称为三阶 Runge-Kutta 方法

48

49

50 因此 ----(8) ----(9) 比较 (7)(8) 两式, 可知 因而三阶 R-K 方法 (7) 具有 3 阶精度 方法 (7)

51 类似于 (7) 式, 还可构造四阶 ( 经典 )Runge=Kutta 方法 (10) 因而方法 (10) 有 4 阶精度

52 例 1. 使用高阶 R-K 方法计算初值问题 解:解: (1) 使用三阶 R-K 方法 (7) (7) RK.m

53 其余结果如下 : (2) 如果使用四阶 R-K 方法 (10) (10) n xn k1 k2 k3 yn

54 其余结果如下 : n xn k1 k2 k3 k4 yn

55

56 四、收敛性与稳定性 一、收敛性 定义１ : 如果一个数值方法对任意固定的点 ， 当 时有 ，则称该法 是收敛的． 定理 1 如果 关于 满足利普希兹条件, 即存在常数 L, 使 (6.27) 且 有界, 则欧拉方法 (6.4) 的整体截断误差满足估计式 (6.28) 其中

57 二、稳定性 绝对稳定域 定义 2: 设用某一数值方法计算 时，所得到的实际结果 为 且由误差 引起以后各节点处 的 误差为 如果总有 则称该法是绝对稳定的． 用试验方程 (6.30) （其中 为常数，可以为复数）并且把能使某一数值方 法绝对稳定的 的允许取值范围称为该法的绝对稳定域。

常微分方程组和高阶 微分方程的数值解法简介 一、常微分方程组的数值解法 下列包含多个一阶常微分方程的初值问题 (1) 称为常微分方程组的初值问题

59 (1) 式具有 n 个未知函数 做如下假设 则 (1) 式化为矩阵形式 (2)

60 只要将以前所介绍的各种求解方法中的函数转化为 函数向量, 即可得到相应的常微分方程组的数值解法 这里只介绍求解微分方程组的计算机实现 [x,Y] = ODE45('F',xspan,Y0) 首先编制微分方程组的函数文件： function z=F(x,Y) z=F(x,Y); 然后使用求解命令 ode45 F 为微分方程组的文件名 xspan 为需求值的节点向量 Y0 为函数向量的初值 x 为自变量， Y 为函数值矩阵

61 例 1. 求解微分方程组 解:解: 首先编制微分方程组的 m 文件 function z=fun(x,y) z(1)=x-y(1)+2*y(2); z(2)=2*x-3*y(1)-5*y(2); z=z'; 再编写求解程序

62 xspan=0:0.1:2; y0=[0,1.5]'; [x,y]=ode45('fun',xspan,y0) plot(x,y) 运行 后得

63 二、高阶常微分方程的数值解法简介 例 2. 求下列高阶微分方程的数值解 解:解: 显然 假设 则

64 即二阶问题化为微分方程组的初值问题 Gaojiefangcheng.m gaojie.m function z=gaojie(x,y) z=[y(2);y(3);y(1)*y(2)+3*y(3)];

65 x y function gaojiefangcheng() xspan=0:0.1:1; y0=[0,1,-1]'; [x,y]=ode45('gaojie',xspan,y0); [x,y(:,1)] plot(x,y(:,1)) xlabel('x') ylabel('y')

66