第九章 常微分方程数值解法

1 第九章常微分方程数值解法 第一节Euler方法 第三节单步法的收敛性和稳定性 第二节Runge Kutta方法 上一页下一页返回 2 上一页下一页返回 本章介绍求解微分方程数值解的基本思想和方法 含有自变量 未知函数和它的一阶导数和高阶导数的方程 常微分方程 它是描述运动 变化规律的重要数学方法之一 分为两类 1 初值问题即给出未知函数及导数在初始点的值 2 边值问题即给出未知函数及 或 它的某些导数在区间两个端点的值 3 考虑一阶常微分方程的初值问题 只要f x y 在 a b R1上连续 且关于y满足Lipschitz条件 即存在与x y无关的常数L使对任意定义在 a b 上的y1 x 和y2 x 都成立 则上述问题解存在唯一解 所谓数值解法就是要计算出初值问题的解函数y x 在一系列离散点a x0 x1 xN b上的近似值 y0 y1 yN 节点间距为步长 通常采用等距节点 即取hi h 常数 yn 称为问题的数值解 数值解所满足的离散方程统称为差分格式 上一页下一页返回 4 第一节欧拉方法 一 欧拉公式 令yn为y xn 的近似值 将上式代入 式可得 此式称为欧拉 Euler 公式 为Euler方法的局部截断误差 上一页下一页返回 5 例1用欧拉公式解初值问题 解 取步长h 0 1 欧拉公式的具体形式为 依次计算可得 y3 y4 y5 y6 y7 y8 y9 y10 上一页下一页返回 6 其部分结果见下表 可见Euler方法的计算结果精度不太高 上一页下一页返回 7 欧拉公式的几何意义 几何意义 用折线近似代替方程的解曲线 因而也称Euler方法为折线法 上一页下一页返回 8 二 后退的欧拉公式 也用一阶差商逼近导数 令yn 1为y xn 1 的近似值 则可得 称为后退Euler公式 已知yn时 必须通过解方程才能求出yn 1 这样的公式称为隐式公式 而Euler公式为显式公式 Euler公式和后退Euler公式都是由yn去计算yn 1 因此 称它们为单步法 上一页下一页返回 9 显然 p越大 精度越高 三 局部截断误差与方法的阶 Euler方法的精度 其中 上一页下一页返回 10 所以 Euler方法具有1阶精度 将 在点 处一阶Taylor展开 上一页下一页返回 11 所以 后退的Euler方法也具有1阶精度 隐式Euler方法的精度 上一页下一页返回 12 显 隐式两种算法的平均 欧拉公式的改进 其局部误差为 此公式具有2阶精度 称平均公式或梯形公式 梯形公式可由下迭代式计算 其中迭代初值是Euler公式提供 上一页下一页返回 13 四 改进的欧拉公式 注 此法亦称为预测 校正法 可以证明该算法具有2阶精度 同时可以看到它是个单步递推格式 比隐式公式的迭代求解过程简单 它的精度高于显式欧拉法 上一页下一页返回 14 为了便于编程 常将改进的欧拉公式写为 上一页下一页返回 15 例2用改进的欧拉法解例1中的初值问题 解 取步长h 0 1 改进欧拉法的具体形式为 具体计算过程如下 上一页下一页返回 16 依次计算可得 y3 y4 y5 y6 y7 y8 y9 y10 其部分结果见下表 上一页下一页返回 17 例3对下面的初值问题 解 1 取步长h 0 1 欧拉方法的具体公式为 2 取步长h 0 1 改进的欧拉方法的具体公式为 取步长h 0 1 分别用Euler方法 改进的Euler方法求数值解 上一页下一页返回 18 计算结果见下表 上一页下一页返回 19 第二节龙格 库塔法 基本思想 考察改进的欧拉法 可以将其改写为 斜率一定取k1k2的平均值吗 步长一定是一个h吗 只要能对平均斜率提供一种近似算法 就能得到一种对应的差分格式 上一页下一页返回 20 例如取m个点的斜率构造如下形式的公式 该公式称为m级龙格 库塔 Runge Kutta 公式 简称R K公式 求解 只需将公式的局部截断误差在xn点进行Taylor展开 令其前面尽可能多的项为0 便可导出ai bij ci所满足的方程组 即可从中求出这些系数 上一页下一页返回 21 以m 2的情形为例说明建立R K公式的方法 其局部截断误差为 上一页下一页返回 22 因此有 而对于h3 若将k2的Taylor展开式多取一项 会发现h3项的系数不可能为0 而对于上式有无穷多个解 它的每一组解都给出了一个局部截断误差为的二级R K公式 即二阶R K公式 这里有个未知数 个方程 3 2 上一页下一页返回 23 常用的标准四阶R K公式 经典R K方法 最常用的四阶标准R K公式 经典R K方法 为 上一页下一页返回 24 例用四阶标准R K公式解初值问题 解 取h 0 2 四阶标准R K法的具体格式如下 上一页下一页返回 25 已知 上一页下一页返回 26 至少具有四位有效数字 比较 上节用改进的Euler公式计算 取h 0 1 最多具有四位有效数字 上一页下一页返回 27 改进的Euler公式每前进一步只要计算两次f值 而4阶R K公式每前进一步要计算四次f值 但改进的Euler法的步长比4阶R K法的小一半 两者计算总量差不多 而4阶R K法的效果要比改进的Euler法好 由于龙格 库塔法的导出基于泰勒展开 故精度主要受解函数的光滑性影响 对于光滑性不太好的解 最好采用低阶算法而将步长h取小 上一页下一页返回 28 第三节单步法的收敛性与稳定性 收敛性 Convergency 例 就初值问题考察欧拉显式格式的收敛性 解 该问题的精确解为 欧拉公式为 对任意固定的x xi ih 有 上一页下一页返回 29 稳定性 Stability 例 考察初值问题在区间 0 0 5 上的解 分别用欧拉法 隐式欧拉法和改进的欧拉格式计算数值解 1 0000 2 00004 0000 8 00001 6000 101 3 2000 101 1 00002 5000 10 16 2500 10 21 5625 10 23 9063 10 39 7656 10 4 1 00002 50006 25001 5626 1013 9063 1019 7656 101 1 00004 9787 10 22 4788 10 31 2341 10 46 1442 10 63 0590 10 7 Whatiswrong 上一页下一页返回 30 一般分析时为简单起见 只考虑试验方程 常数l 0 可以是复数 当步长取为h时 将某算法应用于上式 并假设在初值产生误差 则若此误差以后逐步衰减 就称该算法相对于z lh绝对稳定 z的全体构成绝对稳定区域 我们称算法A比算法B稳定 就是指A的绝对稳定区域比B的大 上一页下一页返回 31 例 考察隐式欧拉法 可见绝对稳定区域为 注 一般来说 隐式欧拉法的绝对稳定性比同阶的显式法的好 上一页下一页返回 32 一些单步法的绝对稳定区间见下表 Euler方法改进的Euler方法三阶R K法四阶R K法隐式Euler法梯形法 2