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

Presentation on theme: "1 第八章常微分方程初值问题的数值解法. 2 第八章 常微分方程数值解法 8.1 引言 ( 基本求解公式 )8.1 引言 ( 基本求解公式 ) 8.2 Runge-Kutta 法8.2 Runge-Kutta 法 8.3 微分方程组和高阶方程解法简介8.3 微分方程组和高阶方程解法简介."— Presentation transcript:

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

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

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

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

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

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

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

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

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

10 10 --------(5) ( 一 ) Euler 公式

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

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

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

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

15 15 得 依此类推, 有 0 1.0000 0.1000 1.1000 0.2000 1.1918 0.3000 1.2774 0.4000 1.3582 0.5000 1.4351 0.6000 1.5090 0.7000 1.5803 0.8000 1.6498 0.9000 1.7178 1.0000 1.7848

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

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

18 18 依此类推, 得 0 1.0000 0.1000 1.0918 0.2000 1.1763 0.3000 1.2546 0.4000 1.3278 0.5000 1.3964 0.6000 1.4609 0.7000 1.5216 0.8000 1.5786 0.9000 1.6321 1.0000 1.6819 比较不同的结果

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

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

21 21

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

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

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

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

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

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

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

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

30 30 (2) 梯形公式

31 31

32 32 (3) 改进 Euler 公式 改进 Euler 公式 x y 0 1.0000 0.1 0.9050 0.2 0.8190 0.3 0.7412 0.4 0.6708 0.5 0.6071 使用 MATLAB 软件 Euler2.m 结果为

33 33 0.9050 0.8190 0.7412 0.6708 0.6071 Euler 公式 梯形公式 改进 Euler 公式 结果比较 Euler 法的精度不如梯形公式

34 34

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

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

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

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

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

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

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

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

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

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

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

46 46 令 令 参照 Simpson 求解公式

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

48 48

49 49

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

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

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

53 53 其余结果如下 : (2) 如果使用四阶 R-K 方法 (10) (10) n xn k1 k2 k3 yn 1.0000 0.1000 1.0000 1.1025 1.2555 1.1111 2.0000 0.2000 1.2345 1.3755 1.5945 1.2499 3.0000 0.3000 1.5624 1.7637 2.0922 1.4284 4.0000 0.4000 2.0404 2.3423 2.8658 1.6664 5.0000 0.5000 2.7768 3.2587 4.1634 1.9993

54 54 其余结果如下 : n xn k1 k2 k3 k4 yn 1.0000 0.1000 1.0000 1.1025 1.1133 1.2351 1.1111 2.0000 0.2000 1.2346 1.3756 1.3921 1.5633 1.2500 3.0000 0.3000 1.5625 1.7639 1.7908 2.0423 1.4286 4.0000 0.4000 2.0408 2.3428 2.3892 2.7805 1.6667 5.0000 0.5000 2.7777 3.2600 3.3476 4.0057 2.0000

55 55

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

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

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

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

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

61 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 62 xspan=0:0.1:2; y0=[0,1.5]'; [x,y]=ode45('fun',xspan,y0) plot(x,y) 运行 后得

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

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

65 65 x y 0 0 0.1000 0.0945 0.2000 0.1754 0.3000 0.2381 0.4000 0.2765 0.5000 0.2822 0.6000 0.2436 0.7000 0.1451 0.8000 -0.0342 0.9000 -0.3230 1.0000 -0.7586 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 66