《常微分方程数值解》PPT课件.ppt

第9章常微分方程数值解法 9 1引言 实际中 很多问题的数学模型都是微分方程 我们可以研究它们的一些性质 但是 只有极少数特殊的方程有解析解 对于绝大部分的微分方程是没有解析解的 常微分方程作为微分方程的基本类型之一 在自然界与工程界有很广泛的应用 很多问题的数学表述都可以归结为常微分方程的定解问题 很多偏微分方程问题 也可以化为常微分方程问题来近似求解 本章讨论常微分方程的数值解法 9 1Euler方法 对于一个常微分方程 通常会有无穷个解 如 因此 我们要加入一个限定条件 通常会在端点出给出 如下面的初值问题 为了使解存在唯一 一般 要加限制条件在f上 要求f对y满足Lipschitz条件 要计算出解函数y x 在一系列节点a x0 x1 xn b处的近似值 节点间距为步长 通常采用等距节点 即取hi h 常数 常微分方程的解是一个函数 但是 计算机没有办法对函数进行运算 因此 常微分方程的数值解并不是求函数的近似 而是求解函数在某些节点的近似值 设解函数在节点的近似为 由差商公式 我们有 则 向前差商公式 可以看到 给出初值 就可以用上式求出所有的 基本步骤如下 解差分方程 求出格点函数 数值方法 主要研究步骤 即如何建立差分方程 并研究差分方程的性质 这种方法 称为数值离散方法 求的是在一系列离散点列上 求未知函数y在这些点上的值的近似 我们的目的 就是求这个格点函数 为了考察数值方法提供的数值解 是否有实用价值 需要知道如下几个结论 步长充分小时 所得到的数值解能否逼近问题的真解 即收敛性问题 误差估计 产生得舍入误差 在以后得各步计算中 是否会无限制扩大 稳定性问题 做等距分割 利用数值微分代替导数项 建立差分方程 1 向前差商公式 所以 可以构造差分方程 称为局部截断误差 显然 这个误差在逐步计算过程中会传播 积累 因此还要估计这种积累 欧拉法 欧拉法的局部截断误差 欧拉法具有1阶精度 欧拉公式的改进 隐式欧拉法 implicitEulermethod 由于未知数yi 1同时出现在等式的两边 不能直接得到 故称为隐式 implicit 欧拉公式 而前者称为显式 explicit 欧拉公式 一般先用显式计算一个初值 再迭代求解 隐式欧拉法的局部截断误差 即隐式欧拉公式具有1阶精度 例5 1用欧拉法求初值问题 当h 0 02时在区间 0 0 10 上的数值解 方程真解 梯形公式 trapezoidformula 显 隐式两种算法的平均 注 有局部截断误差 即梯形公式具有2阶精度 比欧拉方法有了进步 但注意到该公式是隐式公式 计算时不得不用到迭代法 其迭代收敛性与欧拉公式相似 中点欧拉公式 midpointformula 假设 则可以导出即中点公式具有2阶精度 需要2个初值y0和y1来启动递推过程 这样的算法称为双步法 double stepmethod 而前面的三种算法都是单步法 single stepmethod 简单 精度低 稳定性最好 精度低 计算量大 精度提高 计算量大 精度提高 显式 多一个初值 可能影响精度 使用梯形公式计算时 常采用Euler方法提供迭代初值 则梯形法的迭代公式为 分析迭代过程的收敛性 可得 于是有 梯形方法虽然提高了精度 但算法复杂 在应用迭代公式进行实际计算的时候 每迭代一次 都要重新计算函数值 而迭代又要反复进行若干次 计算量很大 而且往往难以预测 为了控制计算量 通常只迭代一两次就转入下一步计算 从而简化算法 改进欧拉法 modifiedEuler smethod 注 此法亦称为预测 校正法 predictor correctormethod 可以证明该算法具有2阶精度 同时可以看到它是个单步递推格式 比隐式公式的迭代求解过程简单 后面将看到 它的稳定性高于显式欧拉法 或者表示成下列平均化形式 上式又称改进Euler公式 也可写成 例1 解该方程为贝努利方程 其精确解为 欧拉公式的具体形式为 改进Euler公式的具体形式为 取步长计算结果见下表 9 2 4 单步法的局部截断误差与阶 定义1 在yn准确的前提下 即yn y xn 时 用数值方法计算yn 1的误差Tn 1 y xn 1 yn 1 称为该数值方法计算时yn 1的局部截断误差 Tn 1之所以称为局部的 是假设在xn前各步没有误差 当yn y xn 时 计算一步 则有 则 Euler法的局部截断误差为 局部截断误差主项 定义2数值方法的局部截断误差为 则称这种数值方法的阶数是P 步长 h 1 越小 P越高 则局部截断误差越小 计算精度越高 局部截断误差主项 9 3龙格 库塔法 Runge KuttaMethod 建立高精度的单步递推格式 单步递推法的基本思想是从 xi yi 点出发 以某一斜率沿直线达到 xi 1 yi 1 点 欧拉法及其各种变形所能达到的最高精度为2阶 斜率一定取K1K2的平均值吗 步长一定是一个h吗 首先希望能确定系数 1 2 p 使得到的算法格式有2阶精度 即在的前提假设下 使得 Step1 将K2在 xi yi 点作Taylor展开 Step2 将K2代入第1式 得到 Step3 将yi 1与y xi 1 在xi点的泰勒展开作比较 要求 则必须有 这里有个未知数 个方程 3 2 存在无穷多个解 所有满足上式的格式统称为2阶龙格 库塔格式 注意到 就是改进的欧拉法 Q 为获得更高的精度 应该如何进一步推广 其中 i i 1 m i i 2 m 和 ij i 2 m j 1 i 1 均为待定系数 确定这些系数的步骤与前面相似 最常用为四级4阶经典龙格 库塔法 ClassicalRunge KuttaMethod 由于龙格 库塔法的导出基于泰勒展开 故精度主要受解函数的光滑性影响 对于光滑性不太好的解 最好采用低阶算法而将步长h取小 例取步长h 0 2 用经典格式求解初值问题 解 由四阶经典龙格 库塔公式可得 9 3 3变步长的龙格 库塔法在微分方程的数值解中 选择适当的步长是非常重要的 单从每一步看 步长越小 截断误差就越小 但随着步长的缩小 在一定的求解区间内所要完成的步数就增加了 这样会引起计算量的增大 并且会引起舍入误差的大量积累与传播 因此微分方程数值解法也有选择步长的问题 以经典的四阶龙格 库塔法为例 从节点xi出发 先以h为步长求出一个近似值 记为 由于局部截断误差为 故有 当h值不大时 式中的系数c可近似地看作为常数 然后将步长折半 即以为步长 从节点xi出发 跨两步到节点xi 1 再求得一个近似值 每跨一步的截断误差是 因此有 这样 由此可得 这表明以作为的近似值 其误差可用步长折半前后两次计算结果的偏差 来判断所选步长是否适当