chap7常微分方程数值解法

7 ODE 1 第7章常微分方程 组 的数值解法 刘东毅天津大学理学院数学系 7 ODE 2 第7章常微分方程 组 的数值解法 主要目标 掌握常微分方程初值问题数值解法的基本理论掌握计算机上的常用算法 主要内容 初值问题计算格式的建立Runge Kutta方法一阶常微分方程组与高阶方程的数值解法 7 ODE 3 第7章常微分方程 组 的数值解法 在科学研究和工程实践中会遇到很多微分方程 虽然从理论上可以证明其解的存在性 但其解的解析表达式往往是很难求解的 或者即使可以写出来 但也难于计算 此时 只能借助数值解来解决问题 常微分方程 组 定解问题是自然科学和工程技术领域中常见的数学模型 本章介绍求解此类问题的基本理论和数值解法 7 ODE 4 定义7 0 1若存在常数L 0 使得对一切的x a b 及y 均有 则称f x y 在D上关于y满足Lipschitz条件 其中L称为Lipschitz常数 7 ODE 5 我们首先给出常微分方程初值问题解的存在惟一性定理 定理假设f x y C D 且关于y满足Lipschitz条件 则一阶常微分方程初值问题 7 0 1 存在唯一解 下面在此前提下 我们讨论上述初值问题 7 0 1 的数值解法 7 ODE 6 然后在节点上建立逼近于原初值问题的计算格式 或差分格式 由此计算出原问题的解y x 在节点x1 x2 xN处的近似值 y1 y2 yN 称它们为常微分方程初值问题的数值解 相邻两个节点的距离hn xn 1 xn称为步长 通常取定步长h 0 即节点xn x0 nh n 0 1 N 其基本思想是在区间 a b 上引入一系列节点 7 ODE 7 7 1初值问题计算格式的建立 1 数值微分方法 在等距节点下讨论问题 利用两点数值微分公式 7 1 1计算格式的建立 将上式代入初值问题 7 0 1 有 7 1 1 7 ODE 8 略去余项 并以数值解yn yn 1代替y xn 及y xn 1 则得差分方程 上式称为Euler公式 利用此式可由初值y0出发按 步进式 方法 逐步求得数值解y1 y2 yN 由于计算yn 1时 只用到它前一步的结果yn 这类公式称为单步法 又因为其关于yn 1是显式形式 故称该Euler公式为显格式 7 ODE 9 如果利用下列数值微分公式 由类似的可导出 上述公式称为后退的Euler公式 此公式为单步法公式 又因为它关于yn 1成隐式形式 所以该公式为隐式公式 简称隐格式 7 ODE 10 类似地 可导出 上述公式称为Euler两步法公式 这因为 当计算yn 1时 要用到yn 1与yn 显然它也是显格式 如果利用下列三点数值微分公式 7 1 3 7 ODE 11 设y x C2 a b 由Taylor公式有 由于故上式即为 略去余项 并以yn yn 1代替y xn 及y xn 1 得到的差分方程正是Euler公式 7 1 4 2 Taylor展开法 7 ODE 12 3 数值积分方法 对 在区间 xn xn 1 上积分 得 则有 对上式中的积分采用不同的数值积分公式可得到不同的差分方程 例如 对上式的积分采用左矩形公式 可得到Euler公式 7 ODE 13 若对此式的积分采用梯形公式 则有 若略去余项 以yn yn 1代替y xn 及y xn 1 得到的差分方程 7 ODE 14 上式称为梯形公式 由于它关于yn 1成隐式形式 故其为隐格式 隐格式求解比较困难 当yn已知时 要求yn 1 需解关于yn 1的非线性方程 在实际应用时 上式常与Euler公式联合使用 构成如下计算格式 7 1 6 7 ODE 15 隐式梯形公式的迭代格式 7 1 7 由上式可以得到一个序列 k 0 1 关于此序列的收敛性 有如下的定理 7 ODE 16 定理7 1 1设f x y 在区域D上关于y满足Lipschitz条件 即 其中L为Lipschitz常数 当步长时 对任意的初值按格式 7 1 7 生成的序列收敛于梯形公式 7 1 6 的解 7 ODE 17 为了减少计算量 可采用预测 校正格式 方法是先用Euler公式求得一个初始近似值称为预测值 再把带入梯形公式右端计算一次求得yn 1称之为校正值 即 上式称为预测 校正公式或改进的Euler公式 上式也可写成如下形式 7 ODE 18 例7 1 1 利用Euler公式与改进的Euler公式求解初值问题 步长h 0 1 解 由步长h 0 1 知节点设数值解为利用Euler公式得 7 ODE 19 计算结果见下表 见书P227表7 1 此初值问题的解析解为 从上表可以看出 数值解yn与解析解y xn 比较 yn精度较差 7 ODE 20 解此问题的改进的Euler公式为 同Euler公式比较 改进的Euler法显然精度提高了 由于误差大小是评价计算格式优劣的重要依据 故需要给出有关误差的概念 计算结果见下表 见书P228表7 2 7 ODE 21 7 1 2截断误差与方法的精度 定义7 1 1称误差en 1 y xn 1 yn 1为数值方法在点xn 1的截断误差 又称整体截断误差 设yk y xk k 0 1 n 则为数值方法在点xn 1的局部截断误差 7 ODE 22 整体截断误差en 1是在没有引进舍入误差的情况下 纯粹因为不准确的计算格式造成的 故又称为方法误差 它不仅与x xn 1这一步的计算有关 而且和xn xn 1 x1这几步的计算都有关系 局部截断误差是假设xn之前各数值解没有误差 仅由xn到xn 1这一步计算由计算格式引起的误差 7 ODE 23 如Euler公式 在点xn 1的整体截断误差en 1 y xn 1 yn 1 局部截断误差 7 ODE 24 定义7 1 2若某数值方法的局部截断误差为则称该方法具有P阶精度 或称其为P阶方法 可以证明 Euler方法的局部截断误差其具有一阶精度 梯形方法的局部截断误差其具有二阶精度 改进的Euler方法的局部截断误差具有二阶精度 7 ODE 25 7 2Runge Kutta方法 继续讨论前面的Taylor展开法 设y x C2 a b 由Taylor公式有 由故上式即为 略去余项 并以yn yn 1代替y xn 及y xn 1 得到Euler公式 7 ODE 26 进一步假设设y x Cp 1 a b 由Taylor公式有 其中 由故式 7 2 1 即为 7 ODE 27 7 2Runge Kutta方法 略去余项 并以数值解yn yn 1代替 7 2 3 中的解析解y xn 及y xn 1 可得到一个差分方程 即 其中余项可写成 注 这里 7 ODE 28 在 7 2 3 中略去余项 用yn yn 1代替y xn 及y xn 1 其中 称 7 2 4 式为求解常微分方程初值问题数值解Taylor的格式 7 ODE 29 由于局部截断误差 可知它是一个p阶方法 当p 1时 上式正是Euler公式 但当p 2时 需要计算f x y x 的高阶导数 特别是对于复杂函数f x y x 的求导 这无疑是大大增加计算量 这是它最大的缺点 因此高阶的Taylor方法是不实用的 德国数学家C Runge及M W Kutta提出了一种改进策略 得到了至今还被作为高精度的单步法广泛使用龙格 库塔法 Runge Kuttamethod 7 ODE 30 7 2 1Runge Kutta方法的基本思想 Runge Kutta方法是利用f在某些点处函数值的线性组合替代 7 2 4 步长h后面括号中的因子来构造差分方程 从而避免了高阶导数的计算 这就是Runge Kutta方法的基本思想 用f在某些点处函数值的线性组合替代这一部分 7 ODE 31 其一般形式为 其中r是上式中调用f的个数 r称为级数 为待定参数 适当确定这些参数 可使之具有尽可能高的精度 如局部截断误差满足 7 ODE 32 7 2 2二阶Runge Kutta方法 考虑r 2的情况 此时有 利用二元函数的一阶Taylor公式 即全微分公式 下面将yn 1与y xn 1 作比较 7 ODE 33 从而有 将上式代入 再由 7 ODE 34 得到 在下面要将yn 1与y xn 1 作比较 使它们的局部截断误差满足 为此考虑y xn 1 7 ODE 35 再根据y xn 1 在点xn的一元3阶Taylor展开式 由刚才已得到的 让它们满足 7 ODE 36 即由 左式含有四个未知元三个方程 因此解不唯一 参数满足左式的一族公式统称二阶Runge Kutta公式 可得参数应满足下列方程组 7 ODE 37 可见 二阶Runge Kutta公式 每计算一步需要两次调用f的函数值 7 ODE 38 7 2 3四阶Runge Kutta方法 当r 4时 类似地可导出四阶Runge Kutta公式 这种公式也有一族 其中常用地有 标准 经典 的Runge Kutta方法 7 ODE 39 Gill公式 Gill公式是标准的Runge Kutta公式的改进形式 这种算法可节省存储单元 并能控制舍入误差的增长 四阶Runge Kutta公式 每一步计算需四次调用f的函数值 计算量较大 但其局部截断误差可达O h5 精度较高 7 ODE 40 例7 2 1用标准的四阶Rung Kutta法解初值问题 取步长h 0 2 解 解此问题的计算公式为 计算结果如下 显然在计算量大致相同的情况下 标准的Runge Kutta方法比改进的Euler方法精确度更高 参见p227和p228的表 7 ODE 41 7 5一阶常微分方程组与高阶方程初值问题的数值解法 7 5 1一阶常微分方程组初值问题 7 ODE 42 写成向量形式 其中 注意 在形式上 7 5 2 与 7 0 1 一样 所以可以把求解常微分方程初值问题的各种数值方法推广到方程组上来 7 ODE 43 利用向量值函数的微积分理论 很容易推导出一阶常微分方程组初值问题的数值解法 如Euler公式 其中 7 ODE 44 7 5 3 的分量形式为 或 7 ODE 45 四阶标准的Runge Kutta公式 设则 7 5 5 的分量形式为 7 ODE 46 四阶标准的Runge Kutta公式的分量形式 其中 7 ODE 47 例7 5 1 试写出用中点公式解下列初值问题的计算公式 解 令 则 再取 由向量形式的中点公式 7 ODE 48 中点公式的向量形式 上述中点公式的分量计算形式为 分量计算形式为 7 ODE 49 整理得分量计算格式 7 ODE 50 7 5 2高阶常微分方程初值问题的数值解法 高阶常微分方程初值问题的一般形式为 引入新变量 则 7 5 7 可化为如下一阶常微分方程组初值问题 7 ODE 51 等价的一阶常微分方程组初值问题 于是 可采用前面7 5 1所介绍的方法求此问题的数值解 7 ODE 52 例7 5 2写出用标准四阶Runge Kutta公式求解 的计算公式 解 令 将此问题转化为一阶常微分方程组初值问题 再令则 7 ODE 53 求此问题数值解的标准四阶Runge Kutta公式为 7