数值分析 第9章 常微分方程数值解法.ppt

第十章 常微分方程数值解法 NumericalMethodsforOrdinaryDifferentialEquations 问题驱动 蝴蝶效应 洛伦兹吸引子 Lorenzattractor 是由MIT大学的气象学家EdwardLorenz在1963年给出的 他给出第一个混沌现象 蝴蝶效应 图10 1 1蝴蝶效应示意图 洛伦兹方程是大气流体动力学模型的一个简化的常微分方程组 该方程组来源于模拟大气对流 该模型除了在天气预报中有显著的应用之外 还可以用于研究空气污染和全球侯变化 洛伦兹借助于这个模型 将大气流体运动的强度x与水平和垂直方向的温度变化y和z联系了起来 参数 称为普兰特数 是规范 化的瑞利数 和几何形状相关 洛伦兹方程是非线性方程组 无法求出解析解 必须使用数值方法求解上述微分方程组 洛伦兹用数值解绘制结果图10 1 1 并发现了混沌现象 1引言 微分方程数值解一般可分为 常微分方程数值解和偏微分方程数值解 自然界与工程技术中的许多现象 其数学表达式可归结为常微分方程 组 的定解问题 一些偏微分方程问题也可以转化为常微分方程问题来 近似 求解 Newton最早采用数学方法研究二体问题 其中需要求解的运动方程就是常微分方程 许多著名的数学家 如Bernoulli 家族 Euler Gauss Lagrange和Laplace等 都遵循历史传统 研究重要的力学问题的数学模型 在这些问题中 许多是常微分方程的求解 作为科学史上的一段佳话 海王星的发现就是通过对常微分方程的近似计算得到的 本章主要介绍常微分方程数值解的若干方法 一 初值问题的数值解法 1 一阶常微分方程初值问题的一般形式 常微分方程的数值解法分为 1 初值问题的数值解法 2 边值问题的数值解法 2 一般构造方法 离散点函数值集合 线性组合结构 近似公式 2 迭代格式的构造 1 构造思想 将连续的微分方程及初值条件离散为线性方程组加以求解 由于离散化的出发点不同 产生出各种不同的数值方法 基本方法有 有限差分法 数值微分 有限体积法 数值积分 有限元法 函数插值 等等 3 如何保证迭代公式的稳定性与收敛性 3 微分方程的数值解法需要解决的主要问题 1 如何将微分方程离散化 并建立求其数值解的迭代公式 2 如何估计迭代公式的局部截断误差与整体误差 称在区域D上对满足Lipschitz条件是指 记 4 相关定义 二 初值问题解的存在唯一性 考虑一阶常微分方程的初值问题 Initial ValueProblem 则上述IVP存在唯一解 只要在上连续 且关于y满足Lipschitz条件 即存在与无关的常数L使 对任意定义在上的都成立 求函数y x 在一系列节点a x0 x1 xn b处的近似值的方法称为微分方程的数值解法 称节点间距为步长 通常采用等距节点 即取hi h 常数 称为微分方程的数值解 三 初值问题的离散化方法 离散化方法的基本特点是依照某一递推公式 值 取 按节点从左至右的顺序依次求出的近似 如果计算需用到前r步的值 则称这类方法为r步方法 2欧拉方法 Euler sMethod 欧拉公式 单步显示公式 向前差商近似导数 亦称为欧拉折线法 Euler spolygonalarcmethod 在假设yi y xi 即第i步计算是精确的前提下 考虑的截断误差Ri y xi 1 yi 1称为局部截断误差 localtruncationerror 定义 若某算法的局部截断误差为O hp 1 则称该算法有p阶精度 定义 欧拉法的局部截断误差 Ri的主项 leadingterm 欧拉法具有1阶精度 例1 用欧拉公式求解初值问题 取步长 解 应用Euler公式于题给初值问题的具体形式为 其中 计算结果列于下表 可用来检验近似解的准确程度 进行计算 数值解已达到了一定的精度 这个初值问题的准确解为 从上表最后一列 我们看到取步长 欧拉公式的改进 隐式欧拉法 implicitEulermethod 向后差商近似导数 由于未知数yi 1同时出现在等式的两边 不能直接得到 故称为隐式 implicit 欧拉公式 而前者称为显式 explicit 欧拉公式 一般先用显式计算一个初值 再迭代求解 隐式欧拉法的局部截断误差 即隐式欧拉公式具有1阶精度 梯形公式 trapezoidformula 显 隐式两种算法的平均 注 的确有局部截断误差 即梯形公式具有2阶精度 比欧拉方法有了进步 但注意到该公式是隐式公式 计算时不得不用到 迭代法 其迭代收敛性与欧拉公式相似 中点欧拉公式 midpointformula 中心差商近似导数 假设 则可以导出即中点公式具有2阶精度 简单 精度低 稳定性最好 精度低 计算量大 精度提高 计算量大 精度提高 显式 多一个初值 可能影响精度 Can tyougivemeaformulawithalltheadvantagesyetwithoutanyofthedisadvantages Doyouthinkitpossible Well callmegreedy OK let smakeitpossible 改进欧拉法 modifiedEuler smethod Step1 先用显式欧拉公式作预测 算出 注 此法亦称为预测 校正法 predictor correctormethod 可以证明该算法具有2阶精度 同时可以看到它 是个单步递推格式 比隐式公式的迭代求解过程 简单 后面将看到 它的稳定性高于显式欧拉法 改进的欧拉法 在实际计算时 可将欧拉法与梯形法则相结合 计算公式为 应用改进欧拉法 如果序列收敛 它的极限便满足方程 改进欧拉法的截断误差 因此 改进欧拉法公式具有2阶精度 例2 用改进Euler公式求解例1中的初值问题 取步长 解 对此初值问题采用改进Euler公式 其具体形式为 计算结果列于下表 改进的Euler法 Euler法 通过计算结果得比较可以看出 改进的Euler方法 的计算精度比Euler方法要高 3龙格 库塔法 Runge KuttaMethod 建立高精度的单步递推格式 单步递推法的基本思想是从 xi yi 点出发 欧拉法及其各种变形所能达到的最高精度为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阶龙格 库塔格式 注意到 就是改进的欧拉法 为获得更高的精度 应该如何进一步推广 其中 i i 1 m i i 2 m 和 ij i 2 m j 1 i 1 均为待定系数 确定这些系数的步骤与前面相似 考虑一阶常微分方程初值问题 将区域 a b 进行分划 若 则 n级显式Runge Kutta方法 二阶Runge Kutta方法 取n 2 记 由此得 另一方面 为使局部截断误差为 应取 改进的Euler方法 取 中点方法 取 二阶Heun方法 取 二级Runge Kutta方法不超过二阶 记 则 因此局部截断误差只能达到 三级Runge Kutta方法 取n 3 记 又由于 因此要使局部截断误差为 必须 Kutta方法 取 三阶Heun方法 取 三级Runge Kutta方法不超过三阶 完全类似于二级Runge Kutta方法的分析 只能达到 三级Runge Kutta方法的局部截断误差 将和都展开到项 易证 四级R K方法 取n 4 局部截断误差为O h5 四阶经典龙格 库塔法 ClassicalRunge KuttaMethod 附注 二阶Runge Kutta方法的局部截断误差只能达到 五阶Runge Kutta方法的局部截断误差只能达到 四阶Runge Kutta方法的局部截断误差只能达到 三阶Runge Kutta方法的局部截断误差只能达到 注 龙格 库塔法的主要运算在于计算的值 即计算的值 Butcher于1965年给出了计算量与可达到的最高精度阶数的关系 由于龙格 库塔法的导出基于泰勒展开 故精 太好的解 最好采用低阶算法而将步长h取小 度主要受解函数的光滑性影响 对于光滑性不 4单步方法的收敛性与稳定性 ConvergencyandStability 收敛性 Convergency 例 就初值问题考察欧拉显式格式的收敛性 解 该问题的精确解为 欧拉公式为 对任意固定的x xi ih 有 稳定性 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 一般分析时为简单起见 只考虑试验方程 testequation 常数 可以是复数 我们称算法A比算法B稳定 就是指A的绝对稳定区