计算方法第6章_常微分方程数值方法

1、第六章 常微分方程数值方法连续问题的离散处理寻求微分方程的解在某些离散点上的值在处的近似值；记号：所求函数在处的（准确）函数值， 计算得到的处的函数（近似）值，节点， 步长； 等步长：；以下若不另作说明，一般总记为等步长。§ 6.1 初值问题的数值方法考察微分方程初值问题：由微分方程理论可知：若函数关于满足条件，即存在与无关的常数，使初值问题的解存在且唯一；6.1.1 法及其变形1、法 由Taylor 展开式： 作局部化假设： ，并略去项，便有 将此式右端作为的近似，便得到公式 Euler公式两者的差（即略去的项）称为局部截断误差，记作：Euler公式也可由其他方法导出，例如由第四章

2、数值导数公式，可有：；解出，并由 替换，便可得Euler公式。又如，根据Newton-leibniz公式：将以“0次插值多项式”替换，即以代入积分，得到数值积分（左矩形）公式，及其误差： 又得到与Taylor 展开式相同的表达式，从而又导出Euler公式。 几何意义折线法若一个公式的局部截断误差为，则称该公式的精度为阶，或该公式为阶公式。Euler公式是1阶公式。注意，以上的截断误差是在局部化假设的前提下得到的，即认定。倘若在每一步都按局部化假设，我们有Euler公式的总体截断误差： 2、后退法若取数值导数公式： 与前相同的推导过程，可以得到 在局部化假设的前提下截去局部截断误差 便得到后退法

3、公式：注意到此公式中的右端也有，需要求解关于的方程才能得到。因此将这类公式称为隐式公式，而将可以通过直接计算得到的公式称为显式公式。后退公式是一阶隐式公式，Euler公式是1阶显式公式6.1.2 多步法1、 梯形公式：在式 右端的积分中，取梯形积分公式，有由此，并据微分方程，可得：梯形公式 局部截断误差： 这是一个2阶隐式公式。2、 Simpson公式：在式 右端的积分中，取Simpson积分公式，有由此，并据微分方程，可得：Simpson公式 局部截断误差： ； 与以前的公式不同，用Simpson公式计算，必须有前2步的函数值：和。因此这种方法称为2步方法，而为启动此算法所需的最初的2个函数

4、值：称为表头。更一般的，若计算必须有至少前2步函数值，则这种方法称为多步法。具体地，若计算必须有前k步的函数值，则这种方法称为k步方法，而为启动该方法所需的最初的k个函数值：称为该方法的表头。与此相对，以前的方法计算，只须前1步的函数值，便称为单步方法。因此，Simpson公式是2步、4阶、隐式方法。3、 Adams方法（线性多步法）在式 右端的积分中，若取具有k+1节点的插值多项式近似替代作为被积函数，导出初值问题的求解方法称为Adams方法。（1）显式Adams方法Adams-Bashforth公式取处的构造插值多项式取代：其中。由于，有由此（以后为方便计，记），可得显式的多步法Adams

5、-Bash forth公式及其局部截断误差 这是步、阶的显式公式。下表是的 的数值（注意：）011123-121223-16532455-5937-947201901-27742616-1274251以下是 的公式推导过程：作变量代换：，以为变量，当的变化区间为时，的变化区间为，且 ，有（2）隐式Adams方法Adams-Moulton公式若取处的构造插值多项式取代，与前一样的方法，可得隐式的多步法Adams-Moulton公式及其局部截断误差 这是步、阶的隐式公式。下表是的 的数值（注意：）：011121121258-1324919-514720251646-264106-19述评：从表可见

6、，对相同的，相同，而，特别是：，而有若干，因而在存在计算误差时，由前步导致的误差显然隐式公式要比显式公式小（显式公式对前步的误差会被放大，而显隐式公式则不会），而且局部截断误差也是隐式公式要比显式公式小，结论：隐式公式的稳定性一般比显式公式好。6.1.3 待定系数法利用 展开式比较有关项的系数，可以直接导出公式待定系数法。例：求以下数值公式的系数使公式具有尽可能高的精度：解：由于，因此由展开式，同时，对所求公式右端各项也作相应的展开，并乘以相应的系数：由于期望尽可能准确，比较各对应项的系数，可得方程组： ，解之可得： ；注意到局部截断误差是，因此 显然，这就是的显式Adams-Bashfort

7、h公式。例：求以下数值公式的系数使公式具有尽可能高的精度：与上例完全一样，比较系数，可得公式： 局部截断误差：这个公式称为Milne公式，稍后我们将用到它。综上所述，从公式的构造过程可见，微分方程初值问题的算法公式基本上源于：Taylor 展开式、数值积分或数值微分，而插值方法一般也是通过数值积分或数值微分完成的。l Taylor 展开式待定系数法直接来自于Taylor 展开式，通过比较系数，形成线性方程组，取得公式与局部截断误差；l 数值积分前面的多步法等，基本上都源于此，例如梯形、Simpson方法等，而Adams方法用的是，当用不同的点生成的不同的插值多项式代换时，便得到不同的公式，包括

8、显式和隐式公式；l 数值微分由 导出Euler公式，而由 可导出后退Euler公式，若取3点数值微分公式可导出更多公式：其中第二个公式就是“中点公式”；当然它们也可以由待定系数法取得：例如第3个公式：用待定系数法求以下隐式公式的系数及局部截断误差：由Taylor 展开式： 比较：为使公式具有尽可能高的精度，比较系数，得方程组：因此，公式是 ，与之相应，有局部截断误差：6.1.4 问题的性态与算法的稳定性1、 问题的性态问题对原始数据的敏感性 原始数据问题的应有的结果原始数据扰动问题的结果问题的条件数 相对误差之比的上确界：由微积分知识：初值问题：固定步长，初值； 计算结果： 因此：由于，当 ；

9、或： 又 得： 因此： 良态 病态例： 方程的解 病态 （图示）例： 方程的解 良态 （图示）例：， 当 病态， 当 良态。（图示）2、 方法的稳定性方法：由初始误差导致的后期误差的可控性。中点法：由数值导数公式 得 可导出中点法 此为2步2阶公式，但稳定性差。例： 真解：取2种不同的表头： （即一步法），下表中列出对不同的步长： 计算的近似值， 按后退法计算获得，按中点法由初值及表头获得，按中点法由初值及表头获得，误差放大则是比值。后退法误差放大5741119091112110798说明：1、 总体误差，由前可知后退法：因此，对于后退法，若，（）可计算得；若，可计算得 显然实际情况良好；而对

10、于中点法，其总体误差 考虑表头都用准确值（注意在浮点数系中运算，仍是有误差的，这可认为是初始误差），若，（）可计算得，而，可计算得，而实际计算并非如此。2、对于同是中点法：由此表可见，两个中点法的不同结果是由不同的表头导致的，实际的不同仅是的不同，而它们之间的误差（或不同）导致的计算结果的误差被放大5000倍甚至1万余倍。说明中点法并不是好的方法，尽管它是一个2阶方法，但实际计算不如后退法这样的1阶方法。定义：步方法，若存在常数其中及是按此方法分别由表头及计算得到，则称此方法是稳定的。此定义因未限制，故实际可要求，因此本定义又“渐近稳定性”或“0稳定性”。实际计算关心对确定步长，一个算法，每步

11、计算误差是否被放大？无法讨论一般方程研究对象：典型方程： 其中：。原因：对维常微分方程组：，为阶矩阵，它的特征值反映了矩阵的性态，当时，方程是良态，若时则是病态（工程中，通常称之为“稳定性”， 时，系统稳定，否则，不稳定）。定义：对于典型方程，及确定的，一个步方法，若由表头及计算得到及，存在常数使则称此方法对是绝对稳定的，全体的集合称为该方法的绝对稳定域。事实上，也可将此定义改写成对于单步法则为： 对典型方程，单步法 总有 例如：法：， 后退法：， 梯形法： ，；由于单步法 因此，当 则方法绝对稳定：法： ，在复平面上 是以 1为中心，1为半径的圆内部；后退法：， 即，在复平面上 是以 1为中

12、心，1为半径的圆的外部。为稍后介绍方便，对一般的步方法，写成即： 例如 法： 后退法： 中点法： 梯形法： 法： ；定义算法的特征多项式：其中 分别称为算法（*）的第一、第二特征多项式。定理：对，若算法的特征多项式的全体零点 都满足，则 使该算法绝对稳定，全体这样的的集合，称为该算法的绝对稳定域。例如：中点法的特征多项式 ，则 的根，由于，可见，中点法总是不稳定的。6.1.5 预估-校正方法原因：显式方法简单；但一般稳定性不好； 隐式方法求解复杂；但一般稳定性好；将两者结合问题：如何组合1、误差估计 两个方法获得同一 的不同近似，以此估计误差（1） 同阶、同步长的两个不同方法 方法a： 方法b

13、：由于步长一般较小，因此可认为 ，（实际上就是估计此的值，因为一旦获得此估计值，便可得到的估计值）将两式相减，得 即 因此 （2） 不同阶、同步长的两个不同方法，设 方法a： 方法b：由于是的高阶小量，所以两者相加（减）时，后者可略去；将以上两式相减，可得 因此： （3） 同阶、不同步长的同一方法，设步长 步长： 步长：仿（1），将两式相减，得 所以 ；2、预估-校正法 （ Predictor-预估, Evaluation-计算, Correct-校正）（1） Heun 方法改进Euler法预估：Euler法 校正：梯形法 保持二阶精度： 证明：记为采用原始的梯形法获得的结果，即：，由本节初所

14、作的假设：函数关于变量满足条件， (2) Adams-Bashforth-Moulton 方法取的Adams-Bashforth显式公式作预估，Adams-Moulton隐式公式校正： 对应局部截断误差公式： ，根据前已介绍的误差估计方法，有 。由于在计算过程中，右端的值均已得到，故在每一步，均可对误差进行监控。在计算中，对于给定的相对误差界，若 此处 Re（Relative error）相对误差界, （Absolute error）绝对误差界， Sm (=Small)是一个与相当的小值，目的是避免时出现的估计失真（实际上是函数的数量级，只须注意）。若相对误差过大，须步长减半（），新的进程中，

15、将由等原有的值计算的近似，而根据本方法，它是的组合，对此中的原未取得的与的值，可采用插值法得到，但应注意，本方法是4阶方法，因此，用插值法取得这些数值时应保证误差也应至少是的，而5点插值多项式的余项： 且 ，因此用此5点插值方法可以实现误差是的。为此，取最近相邻的5点构造插值多项式，可得所求点处的近似值： ；(3) Milne-Simpson方法以Milne公式为预估，Simpson公式做校正。 他们的局部截断误差公式分别为：；与前相同，两者相减有： 可导出预估公式-Milne公式的误差估计 与数值积分的Romberg公式的导出思想一样，现在设法将这部分的误差“归还”给预估的。由于在得到后，希

16、望对进行修正，使修正后的值更接近再进入校正公式参与计算，这样得到的校正值有望更符合Simpson公式的实际效果。但由于尚未计算，只能用 近似替换（实际上，是，即：，因此有 注：这是4步方法，必须有表头才能开始计算。预估步与校正步当便可执行，但在第一次的修正步计算中，由于此前无，因此只能令，也即第一次因为无法修正所以不进行修正。（4） 修正Milne-Hamming方法解初值问题的Hamming公式：以Milne公式为预估，Hamming公式进行校正。与前相同，两者的局部截断误差相减有： 可导出预估公式-Milne公式的误差估计 校正公式-Hamming公式的误差估计： 与前一样，现在设法将误差

17、“归还”给预估的，以及；由于通常将最终结果记作，因此将Hamming校正公式计算结果记作；与前Milne-Simpson方法一样，在修正时，尚未产生，因此只能使用，从而形成公式： 注：仍与Milne-Simpson方法一样，这是一个4步方法，必须有表头才能开始计算。预估步与校正步当便可执行，但在第一次的修正步计算中，由于此前无，因此只能令。注：如果从求解隐式方程的迭代法来看，这校正步可以认为是迭代了一步，而迭代的初值则是预估值，或是预估的一个修正值。因此从迭代收敛的角度看，应对初值提出一定的条件，通过研究，可以得到如下结果，各方法对步长的限制：A-B-M4方法： Milne-Simpson方法

18、： 修正Milne-Hamming方法： 。6.1.6 Runge-Kutta方法将Heun 方法（改进Euler法）改写成如下形式：此为一单步、2阶、显式公式（比较梯形：2阶隐式）。若注意此处的，它们是两个不同的增量，Heun 方法则是从出发加上若干个增量的加权平均值后得到，将此推广到一般形式，出发加上个增量的加权平价取得，便有：-级Runge-Kutta公式为导出该公式，先引入二元函数Taylor展开式：以二级Runge-Kutta公式为例： 先考虑准确值的展开：其中 因此 另一方面，考虑近似值的展开，由二元函数Taylor展开式，二级Runge-Kutta公式中，故应有：将准确值与近似值

19、展开式作比较，应若令 ，则，取，则 ， Heun 方法（改进Euler法）。取，则 变形Euler法：（注意此公式与中点法的区别）以上两公式均为二阶二级Runge-Kutta公式。常用的四阶四级Runge-Kutta公式（称为经典、古典或标准Runge-Kutta公式ClassicalRunge-Kutta）几何意义：若将上述公式改写为： 其中为分别为原式中中相应的函数值。比较在区间 上的数值积分Simpson公式其中 分别表示在区间 的左、中、右端点的值；它们分别有：；因此可简略地得到该公式的误差公式是 误差估计：在实际计算中，由于Runge-Kutta方法是一个高精度的单步方法。因此，通常

20、作为求解常微分方程的一个独立的方法使用。在计算时，常常要对计算值进行误差估计。以四阶四级的标准Runge-Kutta公式为例：若用同一公式的不同步长方法进行误差估计，例如用步长和计算获得的不同近似进行，为此必须计算11个值（其中： 计算4个值，计算8个值，而是重复的）。这比不估计误差增加7个值的计算量。一个比较简单的方法是Runge-Kutta-Fehlberg方法以此作误差估计，在计算过程中仅需计算6个值的计算量。小结：常微分方程求解的理论问题，主要是问题的性态和方法的稳定性。算法包括基本算法的构造（或：形成）及组合。前者就是通过Taylor展开式Runge-Kutta法、待定系数法形成如M

21、ilne法等，或数值积分（例如：Adams方法，梯形法，Simpson法等），或数值微分（或称数值导数）形成；当然同一方法实际上可以通过多种方法导出。后者就是利用前面已得到的基本方法，进行适当组合，一方面利用显式方法简单，作为预估，利用隐式方法比较稳定，作为校正；进一步，同阶的不同方法，或同一方法不同步长的局部截断误差比较，将误差中的估计项 变换成可计算的，得到实际误差估计值，或直接用来比较误差是否满足要求，或用来对对公式作修正。6.1.7 常微分方程组、高阶常微分方程1、 常微分方程组初值问题：常微分方程组初值问题的解法与常微分方程初值问题的解法几乎是一模一样的，唯一区别就在于所有的公式都用

22、向量代替纯量。为此，记向量则可将常微分方程组初值问题改写成简单的向量形式：.对此向量形式的初值问题的解法，只需将以前一维方程的数值方法改为向量形式便可。例如，解常微分方程组的Euler方法，写成向量形式：， 标准Runge-Kutta公式的向量形式：2、高阶方程：对于高阶常微分方程初值问题，一般可以将它改写成常微分方程组，然后用前述的向量方法求解。例如：3阶常微分方程初值问题，令，并记 则原高阶方程就可以改写成如下形式，从而可按标准方法求解，6.2边值问题数值方法简介6.2.1 差分方法讨论常微分方程： ， 其中 ；对于一般方程： 可以通过以下方法转化为上述形式。对方程乘以，得： ，令，使 ，

23、即 从而一般方程等价于在此方程两端同乘以 ，便得 利用的反函数（由知是关于的单调函数，故有反函数），将 转化为 的函数，便可得 ，此即为所讨论的形式。第一边值条件 第二边值条件 第三边值条件 ；1、 差分离散法用数值导数替代导数等分区间 步长： ， 节点：，由 及数值微分：得 略去局部截断误差（Local truncation error）： 注意：；得：在内点 处的离散逼近方程：，或 。求方程的数值解， 的近似尚差2个方程；由边界条件得到：对第一边值条件，直接代入即可。对第三边值条件（因第二边值为其特殊情况）利用数值导数公式：略去截断误差部分项,可得方程: 至此，已有 个方程，可以求解。2、 差分方程求解：对第一边值条件，方程组： 如果将上述系