常微分方程数值解法-1

《常微分方程数值解法-1》由会员分享，可在线阅读，更多相关《常微分方程数值解法-1（65页珍藏版）》请在金锄头文库上搜索。

1、1,四 常微分方程数值解法,2,常微分方程数值解法,引言（常微分方程数值解法概述） 显式欧拉法、隐式欧拉法、二步欧拉法 局部截断误差与精度 改进的欧拉方法 龙格-库塔方法 收敛性与稳定性简述 一阶常微分方程组与高阶常微分方程,3,引言,一阶常微分方程初值问题：,定理：若 f (x, y) 在某闭区域 R ：,上连续，且在 R 域内满足李普希兹 (Lipschitz) 条件，即存在正数 L，使得对于 R 域内的任意两值 y1, y2，下列不等式成立：,则上述初值问题的连续可微的解 y(x) 存在并且唯一。,4,引言（续）,实际生产与科研中，除少数简单情况能获得初值问题的初等解（用初等函数表示的解

2、）外，绝大多数情况下是求不出初等解的。 有些初值问题即便有初等解，也往往由于形式过于复杂而不便处理。 实用的方法是在计算机上进行数值求解：即不直接求 y(x) 的显式解，而是在解所存在的区间上，求得一系列点 xn (n = 0, 1, 2, ) 上解的近似值。,5,欧拉（Euler）方法,方法一 化导数为差商的方法,由于在逐步求解的过程中，y(xn) 的准确值无法求解出来，因此用其近似值代替。 为避免混淆，以下学习简记：,y(xn)：待求函数 y(x) 在 xn 处的精确函数值 yn ：待求函数 y(x) 在 xn 处的近似函数值,6,代入初值问题表达式可得：,根据 y0 可以一步步计算出函数

3、 y = y(x) 在 x1, x2, x3 x4, 上的近似值 y1, y2, y3, y4 , 常微分方程数值解是一组离散的函数值数据，它的精确表达式很难求解得到，但可以进行插值计算后用插值函数逼近 y(x),7,欧拉方法（续）,方法二 数值积分法,同样以近似值 yn 代替精确值 y(xn) 可得：,将微分方程 y = f (x, y) 在区间 xn, xn+1 上积分：,8,欧拉方法的几何意义,x,y,0,9,隐式欧拉法,在数值积分法推导中，积分的近似值取为积分区间宽度与右端点处的函数值乘积，即：,这样便得到了隐式欧拉法：,隐式欧拉法没有显式欧拉法方便,10,二步欧拉法,在数值积分法推导

4、中，积分区间宽度选为两步步长，即积分区间为：xn-1, xn+1，则：,以 y(x) 在 xn -1, xn 上的近似值代替精确值可得：,中矩形公式,11,梯形公式欧拉法,在数值积分法中，如果用梯形公式近似计算 f (x, y) 在区间 xn, xn+1 上的积分，即：,用近似值代替精确值可得梯形公式欧拉法：,上式右端出现了未知项，可见梯形法是隐式欧拉法的一种；实际上，梯形公式欧拉法是显式欧拉法与隐式欧拉法的算术平均。,12,例,用显式欧拉法、隐式欧拉法、梯形法求解初值问题：,取 h = 0.1，计算到 x = 0.5，并与精确解进行比较,解：由已知条件可得：h = 0.1，x0 = 0， y

5、0 = 1， f (x, y) = - y + x + 1,显式欧拉法：,13,例：（续）,隐式欧拉法：,化简得：,梯形公式欧拉法：,14,计算结果：,本题的精确解为：,15,局部截断误差,为了简化分析某常微分方程数值算法的误差，现假设 yn = y(xn)，即在前一步 yn 准确的前提下，估计：,称上述误差 Tn+1 为该常微分方程数值算法的局部截断误差,如果某个常微分方程数值算法的局部截断误差可表示为 O(h p+1)，则称该数值算法的精度是 p 阶,欧拉法的精度为一阶；二步欧拉法的精度为二阶；梯形法的精度为二阶。,16,泰勒展开法,如果初值问题中的 f (x, y) 充分可微，则可将 y

6、(xn+1) 在点 xn 处展开：,如果只保留线性项，忽略 h2 及以上各项，则：,显式欧拉公式,17,局部截断误差的分析,利用泰勒公式展开，比较各算法与展开式的前几项,将 y(xn+1) 在 xn 点处用泰勒公式展开：,显式欧拉法的局部截断误差：,欧拉法,1 阶精度,18,补充：二元函数微分中值定理,19,y(xn+1) 在 xn 点处展开：,隐式欧拉法：,1 阶精度,20,分别将 y(xn+1), y(xn-1) 在 xn 点处用泰勒公式展开：,二步欧拉法的局部截断误差：,二步欧拉法：,2 阶精度,21,梯形公式欧拉法：,y(xn+1) 在 xn 点处展开：,2 阶精度,22,各种欧拉法的

7、比较,23,改进的欧拉法,从上述例子可以看到，梯形法由于具有二阶精度，其局部截断误差比显式欧拉法和隐式欧拉法小，但梯形法实质上是一种隐式算法 显式欧拉法是一个显式算法，虽然计算量较小，但是精度不高 综合两种方法的长处，可以先用显式欧拉法求出 y(xn+1) 的一个粗略近似值，然后用它代入梯形法公式的右端，用梯形法计算 y(xn+1) 的较为精确的近似值。,24,改进的欧拉法（续）,按照上述思想，可以建立如下预报-校正系统：,按以上两式求解常微分方程的算法称为改进的欧拉法，它还可以表示为：,嵌套形式,平均化形式,2 阶精度,25,用改进欧拉法求上例所述的初值问题并与欧拉法和梯形法比较误差的大小。

8、,解：采用改进欧拉法的嵌套形式：,26,计算结果,可见，改进欧拉法的误差数量级与梯形法大致相同，而比欧拉法小得多。,27,改进的欧拉法的意义,改进的欧拉法的平均化形式,y (xn+1) 在点 xn 处的一阶展开式为：,28,改进的欧拉法的几何意义,0,x,y,29,龙格-库塔(Runge-Kutta)方法,改进的欧拉法（2 阶精度）,y (xn+1) 在点 xn 处的一阶泰勒展开式为：,显式欧拉法（1 阶精度）,30,龙格-库塔方法（续）,显式欧拉法用一个点的值 k1 作为 k* 的近似值 改进的欧拉公式用二个点的值 k1 和 k2 的平均值作为 k* 近似值； 改进的欧拉法比显式欧拉法精度高

9、； 在 xn, xn+1 内多预报几个点的 ki 值，并用其加权平均值作为 k* 的近似值从而构造出具有更高精度的计算公式，这就是龙格-库塔方法的基本思想。,31,二阶龙格-库塔方法,以 k1 和 k2 的加权平均来近似取代 k*,为分析局部截断误差，令 yn = y(xn)，由泰勒公式得：,32,补充：二元泰勒展开式,33,用二元泰勒公式展开,将 k1, k2 代入 中可得：,34,二阶龙格-库塔方法（续）,2 阶精度,35,四个未知变量，只有三个方程，有无穷多组解 每组解的构成的龙格-库塔方法均为二阶,二阶龙格-库塔方法即为改进的欧拉方法,变形的欧拉法 中 点 方 法,36,三阶龙格-库塔

10、方法,三阶龙格-库塔方法是用三个值 k1, k2, k3 的加权平均来近似取代 k*,要使三阶龙格-库塔方法具有三阶精度，必须使其局部截断误差为 O(h4) 将 k1, k2, k3 代入 yn+1 的表达式中，在 (xn, yn) 处用二元泰勒公式展开，与 y(xn+1) 在 xn 处的泰勒展开式比较,37,三阶龙格-库塔方法（续）,类似二阶龙格-库塔方法的推导过程，8 个待定系数 c1, c2, c3, a2, a3, b21, b31, b32 应满足：,8 个未知参数，6 个方程，有无穷多组解,库塔公式,38,四阶龙格-库塔方法,类似可以推出四阶龙格-库塔公式，常用的有：,标准四阶龙格

11、-库塔公式,39,四阶龙格-库塔方法（续）,吉尔（Gill）公式,4 阶以上龙格-库塔公式的计算量太大，并且精度不一定提高，有时反而会降低，因此实际应用中一般选用四阶龙格-库塔已足可满足精度要求。,40,用经典四阶龙格-库塔方法求解前例的初值问题，并与改进 欧拉 法、梯形法在 x5 = 0.5 处比较其误差大小,解：采用经典四阶龙格-库塔公式：,41,四阶R-K方法的精度比二阶方法高得多,精确解为：,R-K方法的误差：,改进欧拉法的误差：,梯形法的误差：,42,变步长的龙格-库塔方法,设 y (xn) 在 xn 处的值 yn = y (xn)，当 xn+1 = xn+ h 时 y (xn+1)

12、 的近似值为 ，由于四阶 R-K 方法的精度为 4 阶，故局部截断误差为：,用四阶R-K方法求解初值问题精度较高，但要从理论上给出误差 | y (xn) - yn | 的估计式则比较困难；那么应如何判断计算结果的精度以及如何选择合适的步长 h？ 通常是通过不同步长在计算机上的计算结果进行近似估计。,43,若以 h/2 为步长，从 xn 出发，经过两步计算，得到 y(xn+1) 的近似值,变步长的龙格-库塔方法（续）,以上每步的截断误差约为 cn(h/2)5，于是两步的局部截断误差为：,于是：,整理得：,44,变步长的龙格-库塔方法（续）,记： ，给定的精度要求为 e,D e，反复将步长折半计算

13、，直至 D e，取最终得到的 作为 y(xn+1) 的近似值。,D e，再将步长折半一次计算，最终得到符合精度要求的 y(xn+1) 的近似值。,45,单步法的收敛性,显式单步法可统一写成：,增量函数，仅依赖于函数 f，且仅仅是 xn, yn, h 的函数,求 y = y(x),h 0 时，近似解 是否收敛到精确解,，它应当 是一个固定节点，因 此 h 0 时应同时附 带 n ,46,单步法的收敛性（续）,对于 p 阶的常微分方程数值算法，当 h 0， n 时，是否 yn+1 y(xn+1)？,p 阶算法的局部截断误差为：,显然：,局部截断误差的前提假设是：,局部截断误差 0 并不能保证算法收

14、敛,47,单步法的收敛性（续）,定义：若求解某初值问题的单步数值法，对于固定的 当 h 0 且 n 时，它的近似 解趋向于精确解 y(xn)，即：,则称该单步法是收敛的。,定义：称 y(xn) - yn 为单步法的近似解 yn 的整体截断 误差。,单步法收敛,48,单步法的收敛性（续）,收敛性定理,若某单步法满足以上条件，则该方法是收敛的,则该单步法的整体截断误差为：,若单步法 具有 p 阶精度，且增量函数 关于 y 满足：,Lipschitz 条件：,初值 y0 是准确的,49,假设在前一步 yn 准确的前提下求得的近似值为：,算法精度为 p 阶，局部截断误差：,50,51,即：,若初值是准

15、确的，则 e 0 = 0 ，从而整体截断误差为：,y = e x 为单调增函数，当 时,当 h 0 且 n 时，则,52,单步法的稳定性,在讨论单步法收敛性时一般认为数值方法本身的计算过程是准确的，实际上并非如此：,初始值 y0 有误差 d = y0 - y(x0) 后续的每一步计算均有舍入误差,这些初始和舍入误差在计算过程的传播中是逐步衰减的还是恶性增长就是数值方法的稳定性问题,53,定义：若一种数值方法在节点 xn 处的数值解 yn 的扰动 ，而在以后各节点 ym (m n) 上产生的扰动为 ，如果：,单步法的稳定性（续）,定义：设在节点 xn 处用数值算法得到的理想数值解为 yn，而实际

16、计算得到的近似解为 ，称差值：,为第 n 步的数值解的扰动。,则称该数值方法是稳定的。,54,单步法的稳定性（续）,欧拉法：,由于函数 f (x, y) 的多样性，数值稳定性的分析相当复杂，通常只研究模型方程,考察模型方程：,即：,假设在节点值 yn 上有扰动 n，在节点值 yn+1 上有扰动 n+1，且 n+1 仅由 n 引起（即：计算过程中不再引起新的误差）,55,欧拉法稳定,即：,欧拉法稳定的条件：,针对模型方程： 的显式欧拉法：,化简得：,56,隐式欧拉法：,考察模型方程：,即：,化简为：,假设 yn 上有扰动 ，则 yn+1 的扰动为：,隐式欧拉法稳定,，上式均成立，所以：,隐式欧拉法稳定是恒稳定的,57,一阶常微分