数值分析第九章ppt课件

数值分析 第九章常微分方程的数值解 一 Euler方法 三 单步法的收敛性和稳定性 二 Runge Kutta方法 四 线性多步法 很多科学技术和工程问题常用常微分方程的形式建立数学模型 但是对于绝大多数的微分方程问题 很难或者根本不可能得到它的解析解 本章重点考察一阶方程的初值问题 的数值解法 就是寻求解y x 在一系列离散点 处的近似值的方法 相邻两个节点间的距离称为步长 一 Euler方法 1欧拉公式 由初值条件表示积分曲线从 出发 并在处的切线斜率为 因此可以设想积分曲线在x x0附近可以用切 线近似的代替曲线 切线方程为 当x x1时 代入有 这样得到y x1 的近似值y1的方法 重复上述方法 当x x2时 依次可以计算出x3 x4 处的近似值y3 y4 由此得到Euler公式 由于用折线近似代替方程的解析解 所以Euler方法也称为Euler折线法 例用Euler法计算初值问题的解在x 0 3时的近似值 取步长h 0 1 解 Euler公式的截断误差 局部截断误差 一步Euler公式产生的误差 总体截断误差 Euler公式的累积总误差 欧拉法的局部截断误差 所以欧拉法具有1阶精度 Lipschitiz条件 若存在正数L 使得对一切 x y1 y2有 则称f x y 满足Lipschitiz条件 欧拉法的总体截断误差 那么 设 为局部截断误差 所以 特别当n m 1时 有 总体误差与h是同阶的 上式还说明 当时 有即 也就是说 ym收敛到方程的准确解 后退Euler公式 隐式欧拉法 隐式欧拉公式 利用向后差商近似导数 由于未知数yn 1同时出现在等式的两边 不能直接得到 故称为隐式欧拉公式 而前者称为显式欧拉公式 一般先用显式计算一个初值 再迭代求解 隐式欧拉法的局部截断误差 即隐式欧拉公式具有1阶精度 2梯形公式和改进Euler方法 梯形公式 设y y x 是的解 故 由此得到 用yn来近似y xn 用yn 1来近似y xn 1 得 梯形公式 梯形公式是隐式的 可以用迭代法求解 具有2阶精度 梯形公式的局部截断误差 中点欧拉公式 假设 则可以导出即中点公式具有2阶精度 需要2个初值y0和y1来启动递推过程 这样的算法称为双步法 而前面的三种算法都是单步法 简单 精度低 稳定性最好 精度低 计算量大 精度提高 计算量大 精度提高 显式 多一个初值 可能影响精度 有没有一种方法 既有这些方法的优点 而没有它们的缺点 改进欧拉法 1 先用显式欧拉公式作预测 算出 2 再将代入梯形公式的右边作校正 得到 注 此法亦称为预测 校正法 可以证明该算法具有2阶精度 同时可以看到它是个单步递推格式 比隐式公式的迭代求解过程简单 例用梯形公式求解初值问题 步长h 0 2 解 梯形公式为 于是 整理得 由y 1 y0 2依次可得y1 y2 y3 y4 y5 例用改进欧拉法求解初值问题 要求步长h 0 2 并计算y 1 2 和y 1 4 解 改进欧拉法公式为 即 由y 1 y0 1计算得 二 Runge Kutta方法 建立高精度的单步递推格式 单步递推法的基本思想是从 xn yn 点出发 以某一斜率沿直线达到 xn 1 yn 1 点 欧拉法及其各种变形所能达到的最高精度为2阶 考察改进的欧拉法 可以将其改写为 斜率一定取K1K2的平均值吗 步长一定是一个h吗 首先希望能确定系数 1 2 p 使得到的算法格式有2阶精度 即在的前提假设下 使得 将改进欧拉法推广为 1 将K2在 xn yn 点作Taylor展开 1二阶Runge Kutta方法 2 将K2代入第1式 得到 3 将yn 1与y xn 1 在xn点的泰勒展开作比较 要求 则必须有 所以存在无穷多个解 所有满足上式的统称为2阶Runge Kutta格式 若则 改进的欧拉方法 若则 中点公式 2四阶Runge Kutta方法 其中 i i 1 m i i 2 m 和 ij i 2 m j 1 i 1 均为待定系数 确定这些系数的步骤与前面相似 由于方程的个数少于未知量的个数 所以方程有无穷多个解 可以根据情况得到几种常用的解 即得到相应的四阶公式 最常用为四阶经典龙格 库塔法 也称为标准四阶龙格 库塔公式 Gill公式 2 龙格 库塔法的导出基于泰勒展开 故精度主要受解函数的光滑性影响 对于光滑性不太好的解 最好采用低阶算法而将步长h取小 注 1 龙格 库塔法的主要运算在于计算Ki的值 即计算f的值 Butcher于1965年给出了计算量与可达到的最高精度阶数的关系 例用标准四阶Runge Kutta法求初值问题 在x 0 1处的近似值 取步长为h 0 1 解 所以 那么 例用标准四阶Runge Kutta法求初值问题 在x 0 4处的近似值 取步长为h 0 2 解 所以 而 所以 1单步法的收敛性 三 单步法的收敛性和稳定性 单步法是在计算yn 1时只用到前一步的信息yn 显式单步法的共同特征是它们都是将yn加上某种形式的增量 得出yn 1 计算公式如下 增量函数 Euler方法的增量函数 改进Euler方法的增量函数 则称 为显式单步法在xn 1处的局部截断误差 例 考察欧拉显式格式的收敛性 解 该问题的精确解为 欧拉公式为 对任意固定的x xi ih 有 Tn 1按h展开的第一项 又称为主项 若局部截断误差的展开式写成 则称为局部截断误差的主项 单步法的收敛定理 设单步法具有p阶精度 其增量函数关于y满足Lipschitz条件 即存在常数L 使对任何的及任意的x有 又设初值y0是准确的 即则总体 截断误差是p阶的 也就是 特别的当时 不论n为何值 总有 即方法收敛 在f x y 对y满足Lipschitz条件下 Euler法 改进Euler法和Runge Kutta法的增量函数 都对y满足Lipschitz条件 所以上述结论对这些方法都成立 例设 是求解微分方程的单步法 试求其局部截断误差的主项 并说出它具有几阶精度 解 考虑在xn处的Taylor展式 所以 该方法的局部截断误差的主项是 具有一阶精度 解 考虑在xn处的Taylor展式 所以 该方法的局部截断误差的主项是 具有二阶精度 2单步法的稳定性 收敛性是在假定每一步计算都准确的前提下 讨论步长时 方法的总体截断误差是否趋于零的问题 稳定性是讨论舍入误差的积累能否对计算结果有严重的影响 例 考察初值问题在区间 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 0 0 5 上的解 分别用欧拉显 隐式格式和改进的欧拉格式计算数值解 一般分析时为简单起见 只考虑试验方程 为复数且Re 0 设在节点值yn处有扰动令 那么 于是 反复应用可得 为使 则 可得如下定义 下面讨论已知的几种方法的绝对稳定区间和绝对稳定区域 显式欧拉法 在复平面上的绝对稳定区域是 即以 1为中心 1为半径的圆域 所以 相应的绝对稳定区间是 隐式欧拉法 后退欧拉法 在复平面上的绝对稳定区域是 是以1为中心 1为半径的圆的外域 所以 相应的绝对稳定区间是 即 如果只考虑 0的实数 则相应的绝对稳定区间对于任意的h都成立 所以是无条件稳定的 梯形公式 在复平面上的绝对稳定区域是 也是无条件稳定的 相应的绝对稳定区间是 龙格 库塔法 而显式1 4阶方法的绝对稳定区域为 例设 是求解微分方程的单步法 分析它的稳定性 解 所以绝对稳定区域是 即为复平面的左半平面 在实数域上是无条件稳定的 解 将代入得 即 例讨论求解初值问题的求解 公式 的稳定性 0为实数 所以绝对稳定区域是 所以 因此是条件稳定的 四 线性多步法 在逐步推进的求解过程中 计算yn 1之前已经求出了一系列的近似值y0 y1 yn 如果充分利用前面信息来预测yn 1 则可期望会获得较高的精度 这就是线性多步法的基本思想 1线性多步法的一般公式 最常用的线性多步法公式为 其中为常数 yn k为y xn k 的近似值 fn k f xn k yn k 特别的当时 上式为显式 否则是隐式 若则称该方法具有p阶精度 若则称局部截断误差的主项为为误差常数 例设yn 1 yn 1 2hf xn yn 为求解常微分初值问题的线性二步法 试求该二步公式的局部截断误差主项 和精度 解 由局部截断误差的定义可知 考虑在xn处的Taylor展式 代入可得 所以局部截断误差的主项为 具有二阶精度 解 局部截断误差为 考虑在xn处的Taylor展式 于是 为使方法具有二阶精度则 解得 因此该方法为 局部截断误差的主项为 解 局部截断误差为 考虑在xn处的Taylor展式 所以 解得 局部截断误差的主项为 2Adams外推公式 考虑用r 1个点 xn k f xn k yn k 构造一个r次多项式来近似 的被积函数f x y x 这里用yn k作为y xn k 的近似值 令fn k f xn k yn k 用Newton向后插值公式来构造r次多项式 即 这里 而 将yn 1作为y xn 1 的近似值 因此得到 Adams外推公式 显式格式 当r 0时 为Euler公式 最常用的是r 3情况 Adams外推系数 3Adams内插公式 用xn 1 xn xn r 1为插值节点 构造f x y x 的r次多项式 得到内插公式 隐式格式 最常用的是r 3情况 Adams外推系数 4预报 校正公式 不论单步法或多步法 隐式公式比显式公式稳定性好 但在实际使用隐式公式时 都会遇到两个问题 一个是隐式公式如何能方便地进行计算 另一个是实际计算步长取多大 如隐式梯形公式 每往前推进一步 不必进行多次迭代 而是采用一阶显式Euler