常微分方程初值问题数值解法

常微分方程初值问题数值解法常微分方程初值问题数值解法 常微分方程初值问题常微分方程初值问题 数值解法数值解法 根据给定的初始条件 确定常微分方程惟一解的问题叫常微分方程初值问题 大多数 实际问题难以求得解析解 必须将微分问题离散化 用数值方法求其近似解 一阶常微分方程的初值问题的提法是 求出函数 y x 使满足条件 1 利用数值方法解问题 1 时 通常假定解存在且惟一 解函数 y x 及右端函数 x y 具 有所需的光滑程度 数值解法的基本思想是 先取自变量一系列离散点 把微分问题 1 离散 化 求出离散问题的数值解 并以此作为微分问题解 y x 的近似 例如取步长 h 0 以 h 剖分 区间 b 令 xi ih 把微分方程离散化成一个差分方程 以 y x 表微分方程初值问题 的解 以 yi表差分问题的解 就是近似解的误差 称为全局误差 因此 设 计各种离散化模型 求出近似解 估计误差以及研究数值方法的稳定性和收敛性等构成了数值 解法的基本内容 离散化方法 常用的有三种 基于数值微分的方法 将方程 1 左端的导数用某个一阶数值微 分公式代替 例如在 xn点以 yn 1 yn h 代替 y 即得到欧拉向前公式 2 若在 xn 1点以 yn 1 yn h 代替则得到欧拉向后公式 3 取 2 3 的平均 可导出二阶精度的梯形公式 4 基于泰勒展开的方法 设计一个算法 假定公式中含有某些待定 常数 在函数光滑的假定下 将其按泰勒展开并与微分方程解 y xn h 的展式中 h 的同幂次项相比较 按照给定的精度阶得到待定常数应满 足的一些方程 通过这些方程确定待定常数 即可得到所要的差分公 式 由此法可导出龙格 库塔公式 设计算公式有下列形式 5 i ij为待定常数 取定 N 值 可按上述泰勒展开的方法确定它们 最常用的显式 4 阶龙格 库塔公式为 6 式中 基于函数数值积分的方法 将微分方程的解 y x 代入方程 1 在 子区间 xn ih xn jh 上积分得到公式 7 积分号下是 x 的函数 若用某些结点上 的值的数值积分公式近似这 个积分 便得到各种差分公式 特别地 若取 i 0 j 1 并用 在结点 xn xn 1 xn 2 上的插值代替 7 式中的被积函数 便得亚当斯外推公式 4 阶亚当斯外推公式为 8 若取 在结点 xn 1 xn xn 1 上的插值代替 7 式的被积函数 则得亚当斯 内插公式 4 阶亚当斯内插公式为 9 解法可按计算 yn 1时用多少个结点上的值分为单步法和多步法 又 可以按 yn 1出现的形式分为显式法和隐式法 单步法是指已知结点 xn上 yn的值便可计算 yn 1的值的解法 如 2 3 4 单步法是可以自己起步的 即可从方程的初值 y0一步步算出 y1 y2 的值 多步法是指已知 yn yn 1 yn k 1 k 2 的值才能计算 yn 1的值的解法 又 称 k 步法 例如 8 是四步法 9 是三步法 多步法不能自己起步 即 给了初值 y0以后 还要用其他解法 如单步法 算出 y1 y2 yk 1后 才能 使用多步法 继续往下计算 多步法公式若对 yi和 i都是线性的 则称 作线性多步法 k 步线性多步法的一般形式为 显式法的公式中 未知的 yn 1明显地被表示 即公式中除 yn 1一项 外 其他的项中不再含有 yn 1 如公式 2 隐式法的公式不显含 yn 1 求未知的 yn 1时一般需要解方程 如公 式 3 或 4 通常用各种迭代方法解隐式差分方程 也可采用较简单的预 估 校正方法 如使用梯形公式 4 时 可先用显式公式 2 求得 yn 1的预 估值 代入式 4 的函数 n 1中 再求得 yn 1的值 此法又称改进的欧拉 折线法 数值解法满足相容的 收敛的 数值稳定的条件时 才有实用价 值 为此要研究以下的一些问题 相容性 将微分方程离散化所带来的误差叫截断误差 当 h 0 时 截断误差趋于零 则称离散化后的方程与微分方程具有相容性 表示 离散化后的方程是微分方程的近似 若截断误差的主要项为 Chp 1 则称 截断误差的阶是 p 1 而称该解法是 p 阶的 p 越大表示离散化后的方 程与微分方程近似程度越高 收敛性 是指当 h 0 时 全局误差 i 0 即离散问题的解 yn收 敛于微分问题的解 y x 这是离散解可用的理论基础 p 阶的解法 即 是当 h 0 时 i以 hp的速度收敛 误差估计 对全局误差 i的估计 是应用数值解法时最关心的 问题 先验估计通常只能给出误差的阶 即误差的主要项中步长 h 的 幂次 一般采用事后估计 即在计算的过程中估计误差 例如用理查森 外推法估计误差 外推法也是提高解的精确度的有效方法 数值稳定性 是指计算过程中 某一步上产生的误差一步一步地 传递下去 是衰减 不增或有界 使得传递下来的误差不致于影响数值解 的精度 至少是不会湮没数值解 数值稳定性是常微分方程数值积分 时必须考虑的问题 1956 年 G 达赫尔斯特证明 存在 2k 阶 k 步线性多步法 但数值稳 定的 k 步线性多步法 当 k 为偶数时 其阶不能超过 k 2 当 k 为奇数 时 其阶不能超过 k 1 称为限制性定理 判别一个数值方法的稳定性时 微分方程 10 有较广的代表性 这里 i 0 许多数值稳定性的定 义都以这个方程为基础 通常称它为测验方程 常微分方程初值问题数值解法常微分方程初值问题数值解法 常微分方程初值问题数值解法常微分方程初值问题数值解法 常微分方程初值问题数值解法常微分方程初值问题数值解法 稳定区域 是指将一个数值积分方法应用于测验方程 10 在 h 复平面上使方法数值稳定的区域 欧拉公式 2 的稳定条件 为 1 h 1 其稳定区域是以 2 0 为直径的圆的内部 图 1 龙格 库塔法 6 的稳定区域由条件 11 来确定 实际上所有 4 阶显式龙格 库塔法的稳定区域都是由条件 11 确定 是图 2 图形的内部 欧拉向后公式 3 的稳定区域是以 0 2 为直径的圆的外部 图 3 稳定区域是稳定条件的几何表示 其作用 在于解线性常系数常微分方程组 12 时 若步长 h 取得使所有 h 是 A 的本征值 都落在稳定区域内 用这个步长积分时数值稳定 理论上稳定区域可以用来选择可用的步 长 h 由图 1 图 2 可知 对于 0 的 只要 h 取得足够小 h 就可落在稳定区域内 对于图 3 h 没有限制 稳定区域包含 h 平面 的整个左半平面的方法叫做 A 稳定的 如公式 3 4 是 A 稳定的 刚性方程组 常微分方程组的初值问题为 式中 皆为 n 维向量 假定微分方程右端函数 的雅可比矩阵的特 征值为 且 Re j 0 j 1 2 n 当比值 时 用通常显式公式计算 h j j 1 2 n 不能超过某个量 也就是 h 必须很小 从而大大增加计算时间 这类问题称为刚性问题 因为 A 稳定的方法从数值稳定来说对步长 h 没有限制 适用于刚性方程组 如 欧拉向后公式 3 梯形公式 4 隐式龙格 库塔公式等都是有效的方法 参考书目 P Henrici Discrete variable Methods in Ordinary Differential Equations John Wiley Sons New York 1962 C W 吉尔著 费景高 刘