第六章常微分方程初值问题的数值解法ppt课件.ppt

《第六章常微分方程初值问题的数值解法ppt课件.ppt》由会员分享，可在线阅读，更多相关《第六章常微分方程初值问题的数值解法ppt课件.ppt（98页珍藏版）》请在三一办公上搜索。

1、计 算 方 法,湖南大学电气与信息工程学院,第六章 常微分方程初值 问题的数值解法,计算方法课程组,定理：若 f (x, y) 在某闭区域 R ：,上连续，且在 R 域内满足李普希兹 (Lipschitz) 条件，即存在正数 L，使得对于 R 域内的任意两值 y1, y2，下列不等式成立：,则上述初值问题的连续可微的解 y(x) 存在并且唯一。,一阶常微分方程的初值问题 :,引言,6.1 基本离散方法,6.2 Runge-Kutta方法,6.3 线性多步法,6.4 收敛性与稳定性,6 常微分方程数值解法,6.5 微分方程组的数值处理,6.6 高阶微分方程的数值处理,考虑一阶常微分方程的初值问题

2、 :,例如：,其解析解为：,6.1 基本离散方法,但是, 只有一些特殊类型的微分方程问题能够得到用解析表达式表示的函数解，而大量的微分方程问题很难得到其解析解。 因此，只能依赖于数值方法去获得微分方程的数值解。,例如：,其解析解为：,很难得到其解析解,要计算出解函数 y(x) 在一系列节点 a = x0 x1 xn= b 处的近似值,通常取节点间距 为步长，通常采用等距节点，即取 hi = h (常数)。它适合计算机求解，应用广泛，具有应用价值。, Taylor级数复习,函数 在 点作Taylor级数展开：,这里 x , x0 都可以是任意一点。,Taylor级数复习,则：,6.1.1 Eul

3、er公式,方法一：利用等距分割，数值微分来代替导数项，建立差分格式。,1、向前差商公式（单步法或显式欧拉法）,所以，可以构造差分方程,称为局部截断误差。显然，这个误差在逐步计算过程中会传播，积累。因此还要估计这种由截断误差带来的积累误差,定义:在假设 yi = y(xi)，即第 i 步计算是精确的前提下，考虑的截断误差 Ri = y(xi+1) yi+1 称为局部截断误差。,定义: 若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶精度。,收敛性：,考察局部误差的传播和积累,2、向后差商公式（单步法或隐式欧拉法）,是隐格式，要迭代求解,可以由向前差商公式求出,3、中心差商公式（二步法或

4、中点欧拉公式）,是多步法，具有2次代数精确度，但该迭代格式不稳定,对微分方程,两边积分有:,类似，可以算出其误差估计式：,所以，有,是个隐式的方法，要用迭代法求解具有2次代数精确度,局部截断误差,4、梯形公式（单步法或梯形欧拉法）,14,以近似值 yn 代替精确值 y(xn) 可得显式欧拉公式：,将微分方程 y = f (x, y) 在区间 xn, xn+1 上积分：,方法二： 数值积分法,1、显式欧拉法,欧拉公式有明显的几何意义,依此类推得到一折线,由于每步计算 只用到 ，故也称此为单步法，,或显格式的单步法。,也称欧拉折线法.,就是用这条折线近似地代替曲线,欧拉方法,17,这样便得到了隐式

5、欧拉法：,隐式欧拉法没有显式欧拉法方便,2、隐式欧拉法,将微分方程 y = f (x, y) 在区间 xn, xn+1 上积分：,以 y(x) 在 xn -1, xn 上的近似值代替精确值可得：,将微分方程 y = f (x, y) 在区间 xn-1, xn+1 上积分：,3、二步欧拉法（中心欧拉法）,用近似值代替精确值可得梯形公式欧拉法：,注：上式右端出现了未知项，可见梯形法是隐式欧拉法的一种；实际上，梯形公式欧拉法是显式欧拉法与隐式欧拉法的算术平均。,4、梯形公式欧拉法,将微分方程 y = f (x, y) 在区间 xn, xn+1 上积分：,取 h = 0.1，计算 y(0.5) 并与精

6、确解进行比较,解：由已知条件可得：h = 0.1，x0 = 0， y0 = 1， f (x, y) = - y + x + 1,显式欧拉法：,例题：用显式欧拉法、隐式欧拉法、梯形法求解初值问题：,其精确解为：,隐式欧拉法：,化简得：,梯形公式欧拉法：,计算结果,其精确解为：,称上述误差 Rn+1 为该常微分方程数值算法的局部截断误差,如果某个常微分方程数值算法的局部截断误差可表示为 O(h p+1)，则称该数值算法的精度是 p 阶,欧拉法的精度为一阶；二步欧拉法的精度为二阶；梯形公式欧拉法的精度为二阶。,6.1.2 局部截断误差,为了简化分析某常微分方程数值算法的误差，现假设 yn = y(x

7、n)，即在前一步 yn 准确的前提下，估计,如果只保留线性项，忽略 h2 及以上各项，则：,显式欧拉公式,泰勒展开法,如果初值问题中的 f (x, y) 充分可微，则可将 y(xn+1) 在点 xn 处展开：,将 y(xn+1) 在 xn 点处用泰勒公式展开：,显式欧拉法的局部截断误差：,欧拉法,1 阶精度,局部截断误差的分析,利用泰勒公式展开，比较各算法与展开式的前几项,补充：二元函数微分中值定理,y(xn+1) 在 xn 点处展开：,隐式欧拉法：,1 阶精度,隐式欧拉法的局部截断误差：,分别将 y(xn+1), y(xn-1) 在 xn 点处用泰勒公式展开：,二步欧拉法的局部截断误差：,2

8、 阶精度,梯形公式欧拉法：,y(xn+1) 在 xn 点处展开：,2 阶精度,梯形-截断误差：,6.1.3 欧拉公式的比较与改进：,向后差商近似导数,由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故称为隐式欧拉公式。,上式称为显式欧拉公式；,显式欧拉法,隐式欧拉法,注意：显式和隐式欧拉法误差都具有1次代数精确度,中点欧拉公式（欧拉二步法）,假设 ，则可以导出即中点公式具有 2 阶精度。,需要2个初值 y0和 y1来启动递推过程，这样的算法称为双步法 /* double-step method */，而前面的三种算法都是单步法 /* single-step method */。,简单

9、,精度低,稳定性最好,精度低, 计算量大,精度提高,计算量大,精度提高, 显式,多一个初值, 可能影响精度,32,各种欧拉法的比较,33,从上述例子可以看到，梯形法由于具有二阶精度，其局部截断误差比显式欧拉法和隐式欧拉法小，但梯形法实质上是一种隐式算法显式欧拉法是一个显式算法，虽然计算量较小，但是精度不高综合两种方法的长处，可以先用显式欧拉法求出 y(xn+1) 的一个粗略近似值，然后用它代入梯形法公式的右端，用梯形法计算 y(xn+1) 的较为精确的近似值。,改进欧拉法思想,改进欧拉法(预估校正欧拉法),Step 1: 先用显式欧拉公式作预测，算出,Step 2: 再将 代入隐式梯形公式的右

10、边作校正，得到,注：此法亦称为预测-校正法 /* predictor-corrector method */。可以证明该算法具有 2 阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法。,Euler方法、隐式Euler方法、梯形方法与单步法计算公式的,显式单步法,对应关系,隐式单步法,显式 Euler方法,隐式Euler方法,梯形方法(隐式),6.1.4 总结,算例：分别用Euler公式和改进的Euler公式求解：,取步长 ，计算y(0.5)的近似值,解：欧拉公式：,改进的Euler公式：,算例 分别用显式Euler方法，梯形方法和预估校

11、正Euler方法 初值问题,解：,取 h =0.1，,(1) 显示Euler方法为:,续,算例 分别用显式Euler方法，梯形方法和预估校正Euler方法解初值问题,解：,取 h =0.1，,(2) 梯形方法为：,续,算例 分别用显式Euler方法，梯形方法和预估校正Euler方法解初值问题,解：,取 h =0.1，,显式Euler方法结果为：,(3)预估校正Euler方法：,续,数值例子表明，梯形方法和预估校正Euler方法比显式Euler方法有更好的精度。,续,改进的欧拉法的意义,改进的欧拉法的平均化形式,y (xn+1) 在点 xn 处的一阶展开式为：,改进的欧拉法的几何意义,0,x,y

12、,6.2 龙格 - 库塔(Runge-Kutta)方法,改进的欧拉法（2 阶精度）,y (xn+1) 在点 xn 处的一阶泰勒展开式为：,显式欧拉法（1 阶精度）,龙格-库塔方法（续）,显式欧拉法用一个点的值 k1 作为 k* 的近似值改进的欧拉公式用二个点的值 k1 和 k2 的平均值作为 k* 近似值；改进的欧拉法比显式欧拉法精度高；在 xn, xn+1 内多预报几个点的 ki 值，并用其加权平均值作为 k* 的近似值从而构造出具有更高精度的计算公式，这就是龙格-库塔方法的基本思想。,建立高精度的单步递推格式,6.2.1 二阶龙格-库塔方法,以 k1 和 k2 的加权平均来近似取代 k*,

13、为分析局部截断误差，令 yn = y(xn)，由泰勒公式得：,补充：二元泰勒展开式,用二元泰勒公式展开,将 k1, k2 代入 中可得：,二阶龙格-库塔方法（续）,2 阶精度,四个未知变量，只有三个方程，有无穷多组解由每组解构成的龙格-库塔方法均为二阶,二阶龙格-库塔方法即为改进的欧拉方法,变形的欧拉法中 点 方 法,综上，通过比较对应系数可以得到如下二阶R-K方法：,1、改进的Euler公式,2、变形的欧拉法中点方法,3、Heun公式,6.2.2 三阶龙格-库塔方法,三阶龙格-库塔方法是用三个值 k1, k2, k3 的加权平均来近似取代 k*,要使三阶龙格-库塔方法具有三阶精度，必须使其局

14、部截断误差为 O(h4)将 k1, k2, k3 代入 yn+1 的表达式中，在 (xn, yn) 处用二元泰勒公式展开，与 y(xn+1) 在 xn 处的泰勒展开式比较,三阶龙格-库塔方法（续）,类似二阶龙格-库塔方法的推导过程，8 个待定系数 c1, c2, c3, a2, a3, b21, b31, b32 应满足：,8 个未知参数，6 个方程，有无穷多组解,库塔公式,53,6.2.3 四阶龙格-库塔方法,类似可以推出四阶龙格-库塔公式，常用的有：,标准四阶龙格-库塔公式,54,四阶龙格-库塔方法（续）,吉尔（Gill）公式,4 阶以上龙格-库塔公式的计算量太大，并且精度不一定提高，有时

15、反而会降低，因此实际应用中一般选用四阶龙格-库塔已足可满足精度要求。,一般的RungeKutta法构造,常见的为3阶，4阶公式,例题：用经典四阶龙格-库塔方法求解前例的初值问题，并与改进欧拉法、梯形法在 x5 = 0.5 处比较其误差大小,解：采用经典四阶龙格-库塔公式：,标准四阶龙格-库塔公式,四阶R-K方法的精度比二阶方法高得多,精确解为：,算例： 分别用Euler公式，改进的Euler公式，经典4阶R-K 公式计算一阶常微分方程初值问题。,并与准确解 比较。,解：Euler公式, 改进的Euler公式取步长h=0.1,经典4阶R-K公式取步长h=0.2。,4阶R-K公式：,计算结果见下表

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

17、，于是两步的局部截断误差为：,于是：,整理得：,变步长的龙格-库塔方法（续）,记： ，给定的精度要求为 e,D e，反复将步长折半计算，直至 D e，取最终得到的 作为 y(xn+1) 的近似值。,D e，再将步长折半一次计算，最终得到符合精度要求的 y(xn+1) 的近似值。, 由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h 取小。,6.3 线性多步法,思想: 用若干节点处的 y 及 y 值的线性组合来近似y(xn+1)。,其通式可写为：,当 10 时，为隐式公式; 1=0 则为显式公式。,Adams方法是线性多步法的一个

18、代表，它是利用插值多项式进,行积分得出来的，这样构造线性多步法的方法称为数值求积法，,它是构造线性多步法的一种途径，另外还有Taylor法。,构造线性多步法的途径：,1 数值积分法,显式Adams方法,隐式Adams方法,其它方法,Simpson方法（二步法）,方程(1.1),由Simpson求积公式，得,用 近似 ，得Simpson方法,局部截断误差为,结论：Simpson方法是四阶隐式方法。,例,Milne法,用过点 的插值多项式近似,得 Milne方法,局部截断误差为,结论：Milne方法是四阶显式线性多步法。,例：设,确定式中待定系数0, 1, 2, 0, 1, 2, 3, 使得公式具

19、有4阶精度。,解：,/* y(xi) = yi */,2 Taylor展开法,解：,个未知数个方程,7,5,令 1 = 2 = 0,以 yi+1 取代 yi1，并取 1 = 2 = 0,取 1 = 1, 2 = 0得到辛普森公式,辛普森 公式,例：设,确定式中待定系数0, 1, 2, 0, 1, 2, 3, 使得公式具有4阶精度。,例：确定下列公式：,中的待定系数 ，使公式具有3阶精度。,由于：,只需：,具有3阶精度，但可以进一步验证公式具有4阶精度,例：确定下列公式：,中的待定系数 ，使公式具有3阶精度。,例： 求解,以下用四阶Adams方法来进行求解。,（1）求,（2）求,只要补算,（3）

20、求,只要补算,现列表看用Adams方法求出的误差，精解为,6.4 收敛性与稳定性, 收敛性,例：就初值问题 考察欧拉显式格式的收敛性。,解：该问题的精确解为,欧拉公式为,对任意固定的 x = xi = i h ，有, 稳定性,例：考察初值问题 在区间0, 0.5上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。,1.00002.0000 4.00008.0000 1.6000101 3.2000101,1.00002.5000101 6.25001021.56251023.90631039.7656104,1.00002.50006.25001.56261013.90631019.76

21、56101,1.00004.97871022.47881031.23411046.14421063.0590107,What is wrong ?!,一般分析时为简单起见，只考虑试验方程 /* test equation */,常数，可以是复数,例：考察隐式欧拉法,可见绝对稳定区域为：,注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。,例：隐式龙格-库塔法,而显式 1 4 阶方法的绝对稳定区域为,其中2阶方法 的绝对稳定区域为,无条件稳定,单步法的收敛性,显式单步法可统一写成：,增量函数，仅依赖于函数 f，且仅仅是 xn, yn, h 的函数,h 0 时，近似解是否收敛到精确解,，它应

22、当是一个固定节点，因此 h 0 时应同时附带 n ,单步法的收敛性（续）,对于 p 阶的常微分方程数值算法，当 h 0， n 时，是否 yn+1 y(xn+1)？, 阶算法的局部截断误差为：,显然：,局部截断误差的前提假设是：,局部截断误差 0 并不能保证算法收敛,单步法的收敛性（续）,定义：若求解某初值问题的单步数值法，对于固定的 当 h 0 且 n 时，它的近似 解趋向于精确解 y(xn)，即：,则称该单步法是收敛的。,定义：称 y(xn) - yn 为单步法的近似解 yn 的整体截断 误差。,单步法收敛,单步法的收敛性（续）,收敛性定理,若某单步法满足以上条件，则该方法是收敛的,则该单步

23、法的整体截断误差为：,若单步法 具有 p 阶精度，且增量函数 关于 y 满足：,Lipschitz 条件：,初值 y0 是准确的,假设在前一步 yn 准确的前提下求得的近似值为：,算法精度为 p 阶，局部截断误差：,即：,若初值是准确的，则 e 0 = 0 ，从而整体截断误差为：,y = e x 为单调增函数，当 时,当 h 0 且 n 时，则,单步法的稳定性,在讨论单步法收敛性时一般认为数值方法本身的计算过程是准确的，实际上并非如此：,初始值 y0 有误差 d = y0 - y(x0)后续的每一步计算均有舍入误差,这些初始和舍入误差在计算过程的传播中是逐步衰减的还是恶性增长就是数值方法的稳定

24、性问题,定义：若一种数值方法在节点 xn 处的数值解 yn 的扰动 ，而在以后各节点 ym (m n) 上产生的扰动为 ，如果：,单步法的稳定性（续）,定义：设在节点 xn 处用数值算法得到的理想数值解为 yn，而实际计算得到的近似解为 ，称差值：,为第 n 步的数值解的扰动。,则称该数值方法是稳定的。,单步法的稳定性（续）,欧拉法：,由于函数 f (x, y) 的多样性，数值稳定性的分析相当复杂，通常只研究模型方程,考察模型方程：,即：,假设在节点值 yn 上有扰动 n，在节点值 yn+1 上有扰动 n+1，且 n+1 仅由 n 引起（即：计算过程中不再引起新的误差）,欧拉法稳定,即：,欧拉

25、法稳定的条件：,针对模型方程：的显式欧拉法：,化简得：,隐式欧拉法：,考察模型方程：,即：,化简为：,假设 yn 上有扰动 ，则 yn+1 的扰动为：,隐式欧拉法稳定,，上式均成立，所以：,隐式欧拉法稳定是恒稳定的,6.5 一阶常微分方程组,显式欧拉法,隐式欧拉法,梯形公式欧拉法,以m=2为例,改进的欧拉公式,四阶龙格-库塔公式,6.6 高阶微分方程的初值问题,一般通过引入新的变量，将高阶微分方程化为一阶微分方程组的方法进行求解,m 阶常微分方程：,用变量替换可以将上述的高阶微分方程转化为一阶微分方程组,设：,初始条件为：,则 m 阶微分方程转化为如下的一阶微分方程组：,转化例子,将 转化为微分方程组。,解：令 ，有：,初始条件：,