1、第5章常微分方程数值解法 5 1引言包含自变量 未知函数及未知函数的导数或微分的方程称为微分方程 在微分方程中 自变量的个数只有一个 称为常微分方程 自变量的个数为两个或两个以上的微分方程叫偏微分方程 微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数 如果未知函数y及其各阶导数都是一次的 则称它是线性的 否则称为非线性的 在 常微分方程 中 对于常微分方程的求解 给出了一些典型方程求解析解的基本方法 如可分离变量法 常系数齐次线性方程的解法 常系数非齐次线性方程的解法等 但能求解的常微分方程仍然是有限的 大多数的常微分方程是不可能给出解析解 譬如 这个一阶微分方程就不能用初等函数及其
2、积分来表达它的解 再如 方程 的解 虽然有表可查 但对于表上没有给出的值 仍需插值方法来计算 从实际问题当中归纳出来的微分方程 通常主要依靠数值解法来解决 本章主要讨论一阶常微分方程初值问题 5 1 在区间a x b上的数值解法 数值方法的基本思想对常微分方程初值问题 5 1 式的数值解法 就是要算出精确解y x 在区间 a b 上的一系列离散节点处的函数值的近似值 相邻两个节点的间距称为步长 步长可以相等 也可以不等 本章总是假定h为定数 称为定步长 这时节点可表示为 数值解法需要把连续性的问题加以离散化 从而求出离散节点的数值解 对常微分方程数值解法的基本出发点就是离散化 其数值解法有个基
3、本特点 它们都采用 步进式 即求解过程顺着节点排列的次序一步一步地向前推进 描述这类算法 要求给出用已知信息计算的递推公式 建立这类递推公式的基本方法是在这些节点上用数值积分 数值微分 泰勒展开等离散化方法 对初值问题中的导数进行不同的离散化处理 差分格式 对于初值问题的数值解法 首先要解决的问题就是如何对微分方程进行离散化 建立求数值解的递推公式 递推公式通常有两类 一类是计算yn 1时只用到xn 1 xn和yn 即前一步的值 因此有了初值以后就可以逐步往下计算 此类方法称为单步法 另一类是计算yn 1时 除用到xn 1 xn和yn以外 还要用到即前面k步的值 此类方法称为多步法 5 2Eu
4、ler方法5 2 1Euler格式Euler方法是解初值问题的最简单的数值方法 初值问题的解y y x 称为它的积分曲线 积分曲线上每一点处的切线的斜率等于函数在这点的值 即有 在点xn处为 取 则有 设y xn 的近似值yn已知 用它带入上式右端进行计算 并取计算结果yn 1作为y xn 1 的近似值 则可得 5 2 Euler格式 Euler法的求解过程是 从初始点P0 即点 x0 y0 出发 作积分曲线y y x 在P0点上切线 其斜率为 与x x1直线 相交于P1点 即点 x1 y1 得到y1作为y x1 的近似值 如上图所示 过点 x0 y0 以f x0 y0 为斜率的切线方程为当时
5、 得 这样就获得了P1点的坐标 同样 过点P1 x1 y1 以斜率作直线交x x2于P2点 线段所在的直线方程为 当时 得 由此获得了P2的坐标 当时 得 重复以上过程 就可获得一系列的点 P1 P1 Pn 对已求得点以为斜率作直线 取 从图形上看 就获得了一条近似于曲线y y x 的折线 这样 从x0逐个算出对应的数值解 Euler格式算法实现 1 计算步骤 输入并计算出 使用Euler格式进行计算 输出 并使转到 直至n N结束 2 Euler格式的流程图 例5 1用Euler格式解初值问题 取步长h 0 1 Euler格式的具体形式为 计算公式的精度 常以Taylor展开为工具来分析计算
6、公式的精度 为简化分析 假定yn是准确的 即在的前提下估计误差 局部截断误差 Euler格式的局部截断误差 由 从而 对 在点xn 1处为 取 则有 设y xn 的近似值yn已知 用它带入上式右端进行计算 并取计算结果yn 1作为y xn 1 的近似值 则可得 5 3 后退的Euler格式 5 2 2后退的Euler格式 隐式方程 5 3 常用迭代法求解 用它带入 5 3 的右端 使之转化为显式 直接计算得 然后再用带入 5 3 的右端 又有 如此反复进行迭代 得 后退的Euler格式的局部截断误差 假设 在xn点Taylor展开 从而有 移项整理得 可得 后退的Euler格式的局部截断误差
7、Euler格式的局部截断误差 5 2 3梯形格式比较Euler格式和后退的Euler格式的误差公式 可以看到 如果将这两种方法进行算术平均 即可消除误差的主要部分而获得更高的精度 即 5 4 梯形方法 梯形方法的局部截断误差 假设 在xn点Taylor展开 从而有 移项整理得 可得 梯形方法 5 4 是隐式的 可用迭代法求解 梯形法的迭代公式为 注 迭代过程是收敛的 5 5 5 2 4改进的欧拉公式显式Euler公式计算工作量小 但精度低 梯形公式虽提高了精度 但为隐式公式 需用迭代法求解 计算工作量大 综合Euler公式和梯形公式便可得到改进的Euler公式 先用欧拉公式 5 2 求出一个初
8、步的近似值称为预测值 它的精度不高 再用梯形公式 5 4 对它校正一次 即迭代一次 求得yn 1 称为校正值 这种预测 校正方法称为改进的Euler格式 预测 校正 它可以表示为嵌套形式 或者表示成下列平均化形式 5 6 改进Euler格式算法实现 1 计算步骤 输入并计算出 使用改进Euler格式进行计算 输出 并使转到 直至n N结束 2 改进Euler格式的流程图 例5 2用改进的Euler方法解初值问题 取步长h 0 1 解 改进Euler方法的具体形式 5 2 5Euler两步格式改进的Euler格式中 下面构造与梯形方法在精度上相匹配的显式方法 在点xn处为 则有 设y xn 1
9、y xn 的近似值yn 1 yn已知 用它带入上式右端进行计算 并取计算结果yn 1作为y xn 1 的近似值 则可得 5 7 Euler两步格式 取 局部截断误差 设 前面介绍过的数值方法 无论是Euler方法 后退的Euler方法 梯形方法 还是改进的Euler方法 它们都是单步法 其特点是在计算yn 1时只用到前一步的信息yn可是公式 5 7 中除了yn外 还用到更前一步的信息yn 1 即调用了前两步的信息 故称其为Euler两步法 5 8 由前面知对于 5 8 5 9 对于 5 9 假定则 从而有 可得事后估计式 可以期望 利用这样估计出的误差作为计算结果的一种补偿 有可能使精度得到改
10、善 设以pn和cn分别表示第n步的预估值和校正值 预测 改进 校正 改进 设以pn和cn分别表示第n步的预估值和校正值 预测 改进 校正 改进 计算 计算 注 算法开始之前必须先给出y0 y1和p1 c1 5 3Runge Kutta方法 定义5 1 如果一种方法的局部截断误差为则称该方法具有p阶精度 5 3 2Runge Kutta方法的基本思想 则 由微分中值定理有 y x f x y x 设 称K 为区间 xn xn 1 上的平均斜率 注 只要对K 提供一种算法 就相应的导出一种计算格式 Euler格式可改写成 改进的Euler格式 可改写成 这个处理过程提示我们 如果设法在 xn xn
11、 1 内多预测几个点的斜率值 然后将它们加权平均作为平均斜率K 则有可能构造出具有更高精度的计算格式 Runge Kutta方法的基本思想 5 3 3二阶Runge Kutta方法在上取两点xn和以这两点处的斜率值K1和K2的加权平均 或称为线性组合 来求取平均斜率K 的近似值K 即 式中 K1为xn点处的切线斜率值 K2为点处的切线斜率值 比照改进的Euler法 将视为 即可得 这样设计出的计算格式具有如下形式 上式含有三个待定系数下面通过适当选取这些系数的值 使得格式具有二阶精度 5 10 将以上结果代入 得 设则 其中 在 xn yn Taylor展开 比较系数后可知 只要 成立 格式
12、5 10 的局部截断误差就等于 二阶Runge Kutta格式 5 11 1 若取则这是无穷多解中的一个解 将以上所解的值代入式 5 10 并改写可得 2 若取则此时二阶Runge Kutta法的计算公式为 此计算公式称为变形的Euler格式 5 3 4三阶Runge Kutta方法 为了进一步提高精度 设除外再增加一点 并用三个点的斜率K1 K2 K3加权平均得出平均斜率K 的近似值 这时计算格式具有形式 为了预报点的斜率值K3 在区间内有两个斜率值K1和K2可以用 可将K1 K2加权平均得出上的平均斜率 从而得到的预报值 这样设计出的计算格式具有如下形式 可以通过适当选取这些系数的值 使得
13、格式具有三阶精度 5 12 从而 和前面的推导类似 只要 成立 格式 5 12 的局部截断误差就等于 三阶Runge Kutta格式 5 13 下列是其中的一种 称为Kutta公式 5 3 5四阶Runge Kutta方法 如果需要再提高精度 用类似上述的处理方法 只需在区间上用四个点处的斜率加权平均作为平均斜率K 的近似值 构成一系列四阶Runge Kutta塔公式 具有四阶精度 即局部截断误差是 龙格 库塔方法的推导基于Taylor展开方法 因而它要求所求的解具有较好的光滑性 如果解的光滑性差 那么 使用四阶龙格 库塔方法求得的数值解 其精度可能反而不如改进的欧拉方法 在实际计算时 应当针
14、对问题的具体特点选择合适的算法 9 3 5四阶龙格 库塔法算法实现 1 计算步骤 输入 h N 使用龙格 库塔公式 9 20 计算出y1 输出 并使转到 直至n N结束 2 四阶龙格 库塔算法流程图 程序实现 四阶龙格 库塔法计算常微分方程初值问题 例5 4取步长h 0 2 用经典格式求解初值问题 解 由四阶龙格 库塔公式可得 可同样进行其余yi的计算 本例方程的解为 从表中看到所求的数值解具有4位有效数字 5 3 6变步长的龙格 库塔法在微分方程的数值解中 选择适当的步长是非常重要的 单从每一步看 步长越小 截断误差就越小 但随着步长的缩小 在一定的求解区间内所要完成的步数就增加了 这样会引
15、起计算量的增大 并且会引起舍入误差的大量积累与传播 因此微分方程数值解法也有选择步长的问题 以经典的四阶龙格 库塔法 5 14 为例 从节点xn出发 先以h为步长求出一个近似值 记为 由于局部截断误差为 故有 当h值不大时 式中的系数c可近似地看作为常数 然后将步长折半 即以为步长 从节点xn出发 跨两步到节点xn 1 再求得一个近似值 每跨一步的截断误差是 因此有 这样 由此可得 这表明以作为的近似值 其误差可用步长折半前后两次计算结果的偏差 来判断所选步长是否适当 当要求的数值精度为 时 1 如果 反复将步长折半进行计算 直至 为止 并以上一次步长的计算结果作为 这种通过步长加倍或折半来处
16、理步长的方法称为变步长法 表面上看 为了选择步长 每一步都要反复判断 增加了计算工作量 但在方程的解y x 变化剧烈的情况下 总的计算工作量得到减少 结果还是合算的 Euler方法 后退的Euler方法 改进的Euler方法 显式单步法 隐式单步法 显式单步法 一阶常微分方程初值问题 O h2 O h2 O h3 二阶R K方法 三阶R K方法 显式单步法 显式单步法 Runge Kutta方法 O h3 O h4 经典的四阶R K方法 显式单步法 O h5 5 4单步法的收敛性及稳定性 5 4 1单步法的收敛性 常微分方程 差分方程 是否合理 数值解法的基本思想 定义5 2 若一种数值方法对
17、于任意固定的xn x0 nh 有则称该方法是收敛的 例5 5 就简单的初值问题考查Euler方法的收敛性 求解方程的Euler格式为 从而 解 该问题的准确解为 得 由 收敛 下面考查一般的单步法 显式单步法的计算公式均可写为 其中称为增量函数 注 不同的显式单步法对应不同的增量函数 Euler格式 改进的Euler格式 5 15 单步法的收敛性定理 定理5 1 假设单步法 5 15 满足以下条件 2 增量函数关于y满足Lipschitz条件 1 具有p阶精度 即局部截断误差为O hp 1 3 初值y0是准确的 即y0 y x0 则其整体截断误差 注 从而判断单步法的收敛性 归结为验证增量函数
18、是否满足Lipschitz条件 1 对于Euler方法 由于故当满足Lipschitz条件时它是收敛的 2 对于改进的Euler方法 假定f x y 关于y满足Lipschitz条件 记其Lipschitz常数为L 则 设限定上式表明关于y满足Lipschitz条件 其Lipschitz常数为 收敛 3 类似可以验证Runge Kutta方法的收敛性 练习 验证二阶Runge Kutta方法的收敛性 设限定上式表明关于y满足Lipschitz条件 其Lipschitz常数为 收敛 5 4 2单步法的稳定性稳定性在微分方程的数值解法中是一个非常重要的问题 因为微分方程初值问题的数值方法是用差分格
19、式进行计算的 而在差分方程的求解过程中 存在着各种计算误差 这些计算误差如舍入误差等引起的扰动 在传播过程中 可能会大量积累 对计算结果的准确性将产生影响 这就涉及到算法稳定性问题 例5 6 考察初值问题在区间 0 0 5 上的解 分别用Euler方法 后退的Euler方法和改进的欧拉方法计算数值解 1 0000 2 00004 0000 8 000016 0000 32 0000 1 00002 5000 10 16 2500 10 21 5625 10 23 9063 10 39 7656 10 4 1 00002 50006 250015 62539 062597 6563 1 0000
20、4 9787 10 22 4788 10 31 2341 10 46 1442 10 63 0590 10 7 What swrong 定义5 3 若一种数值方法在节点值yn上产生大小为 的扰动 如果在其后的各节点值ym m n 上产生的偏差均不超过 则称该方法是稳定的 稳定性不仅与算法有关 而且与方程中函数f x y 有关 讨论起来比较复杂 为简单起见 通常只针对模型方程 为保证微分方程本身的稳定性 假定 一般方程若局部线性化 也可化为上述形式 模型方程相对比较简单 若一个数值方法对模型方程是稳定的 并不能保证该方法对任何方程都稳定 但若某方法对模型方程都不稳定 也就很难用于其他方程的求解
21、可见 要保证误差是不增长的 只要选取h充分小 使得 Euler方法的稳定性 模型方程的Euler公式为 设在节点值yn上有一扰动值它的传播使节点值yn 1产生大小为的扰动值 假设用按照Euler格式得出的计算过程不再有新的误差 则扰动值满足 条件稳定 假设用按照后退的Euler格式得出的计算过程不再有新的误差 则扰动值满足 后退的Euler方法的稳定性 模型方程的后退的Euler公式为 设在节点值yn上有一扰动值它的传播使节点值yn 1产生大小为的扰动值 恒稳定 无条件稳定 则 由于 从而恒成立 例5 7 考察初值问题其准确解y e 100 x是一个按指数曲线衰减得很快的函数 用Euler方法
22、解方程 得 1 若取h 0 025 则Euler方法的具体形式为 2 若取h 0 005 则Euler方法的具体形式为 3 若取h 0 025 则后退的Euler方法的具体形式为 用后退的Euler方法解方程 得 解 将三种情况的计算结果列表如下 0 0250 0500 0750 100 精确解 改进Euler格式h 0 025 Euler格式h 0 005 Euler格式h 0 025 节点 1 52 25 3 3755 0625 0 50 250 1250 0625 0 28578 1633 10 22 3324 10 26 6639 10 3 8 2085 10 26 7380 10 35 5308 10 44 5399 10 5 不稳定 稳定 稳定 注 稳定性不但与方法有关 也与步长h的大小有关 当然也与方程中的f x y 有关 例5 8 考察初值问题分别用Euler方法和后退的Euler方法求解 从稳定性角度出发 h应如何取值 解 Euler方法的稳定条件为 则若取h 0 025 计算过程不稳定 若取h 0 005 计算过程稳定 后退的Euler方法是恒稳定的 因此不论h取何值 计算过程都是稳定的
