第八章基于MATLAB的科学计算—常微分方程数值解法

1、科学计算理论、方法及其基于MATLAB的程序实现与分析微分方程（组）数值解法1 常微分方程初值问题的数值解法微分方程（组）是科学研究和工程应用中最常用的数学模型之一。如揭示质点运动规律的Newton第二定律：（）和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题：（）的解为：（）但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：、绝大多数的常微分方程和常微分方程组得不到（有限形式的）解析解；、实际应用中往往只需要知道常微分方

2、程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）；如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。一般的一阶常微分方程（组）的初值问题是指如下的一阶常微分方程（组）的定解问题：（）其中（） （）常微分方程（组）的初值问题通常是对一动态过程（动态系统、动力系统）演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。.1 常微分方程（组）的Cauch问题数值解法概论假设要求在点（时刻），处初值问题

3、（）的解的（近似）值,如果已求得时刻的值或它的近似值(如时刻的值)，那么将式（）的两端在区间上积分(10)可得 (11)或(12)显然,为了利用式（11）或（12）求得的精确值（近似值），必须计算右端的积分，这是问题的关键也是难点所在，如前所述，一般得不到精确的公式解，因此需要采用数值积分的方法求其近似解,可以说,不同的式值积分方法将给出不同的Cauch问题的数值解法。1.2最简单的数值解法Euler 方法假设要求在点（时刻），处初值问题（）的解的近似值。首先对式（）的两端积分，得（13）对于式（13）的右边，如果用积分下限处的函数值代替被积函数作积分（从几何上的角度看，是用矩形面积代替曲边梯

4、形面积），则有（14）进而得到下式给出的递推算法Euler 方法（15）例用Euler 方法解如下初值问题，取，解：由（15）得结果如下：open Euler_Method.m 如果取，其结果如下图所示：Euler_Method 1.3改进的Euler 方法 对于（15）的右边，如果被积函数用积分限和处的函数值的算术平均值代替（几何上，是用梯形面积代替曲边梯形面积），则有（16）进而得到下式给出的递推算法：（17）通常算法（17）比Euler 方法（15）的精度高，但是，按算法（17）求时要解（非线性）方程（组），这是算法（17）不如Euler 方法的方面，为了） 尽可能地保持算法（17）精度

5、高的优点；） 尽可能地利用Euler 方法计算简单的长处；人们采取了如下的称之为改进的Euler 方法的折衷方案：预测 （18）修正 （19） 例Euler 方法与改进的Euler 方法的比较下图是当时比较的结果：open Improved_Euler_Method.m 1.4Euler 方法和改进的Euler 方法的误差分析由Taylor 公式 （19）说明Euler 方法的截断误差是，类似地，由 （20） （21）以及 (22)让式（20）的两端减式（21）的两端，可得 （23）从上述推导Euler 方法、改进的Euler 方法的过程以及例、例容易看出，改进的Euler 方法Euler 方

6、法的精度高，其原因在于：1 在推导Euler 方法时，我们是用待求解函数在一点处的变化率代替在区间上的平均变化率： （24）2 而在推导改进的Euler 方法时，我们是用待求解函数在两点处变化率的平均值代替在区间上的平均变化率；显然，通常比更接近于在区间上的平均变化率。由此启发人们：适当地选取区间上函数若干点处的变化率，用它们加权平均值代替在区间上的平均变化率，近似解的精度应更高。下面将要介绍的RungeKutta法就是基于上述想法得到的。2RungeKutta法RungeKutta法是按选取区间上函数变化率的个数的多少和截断误差的阶数来区分的一系列方法，如 二阶的RungeKutta法（改进

7、的Euler 方法） （25） 三阶的RungeKutta法 （26） 四阶的RungeKutta法 1) 古典形式 （27） 2) Gill 公式(具有减小舍入误差的优点) （28）4 RungeKutta法的一般形式 (29)其中称为增量函数(Increment Function)以及 （30）需要特别指出的是：在确定（29）、（30）中的参数，和时，应该使（29）的右端和适当阶的（如阶）Taylor展式一致，这样，至少对于底阶的R-K法来说，参与加权的斜率个数与方法的阶数是一致的。例如：一阶的R-K法即Euler法，局部的截断误差（Truncation Error）是，所以，当微分方程的

8、阶是一次函数时是精确的；二阶R-K法即改进的Euler法，局部的截断误差是，所以，当微分方程的阶是二次函数时是精确的；类似地，四阶R-K法局部的截断误差是。但是，一般五阶以及五阶以上的R-K法的局部截断误差不再具有上述的特点，因此用“性价比”来衡量，四阶R-K得到了广泛的应用。下面是Butchers Fifth-Order R-K 方法（1964）： （28） （29）4 线性多步法与常微分方程数值解法的分类 一般的常微分方程初值问题的数值解法都是以递推的形式给出的，即是递推算法，根据递推算法，在计算时，已经得到了前面各时刻的近似值，前面介绍的各种数值解法都有一个共同点：在计算时，只用到了前一

9、个时刻（当前时刻）的“信息”预测未来某时刻系统的“状态”，这样的数值解法称为单步法，对于一个动态过程在时刻的状态而言，不仅前一个时刻的信息对它有影响，前若干个时刻的信息通常对它也有影响，显然，单步法的缺点是没有充分利用已得到的信息。基于上述考虑，适当取前若干个时刻的信息，并用的线性组合代替在区间上的平均变化率所给出的方法称为线性多步法。4.1 线性多步法: 从二步法谈起 在改进的Euler方法中预测 （48）较正 （49） 如果改进预测的手段，如利用数值微分的研究成果： (50)那么，由于（48）的局部截断误差是，精度高于Euler公式的局部截断误差，所以，下列预测-校正方法Predictor

10、: （51）Corrector: （52） 尽管步长增大了，精度却提高了。公式（49）、（50）给出的数值方法称为非自始的（non-self-starting）二步法，其原因在于：利用该方法求解微分方程初值问题时，自身不能提供所需的初值。4.3 线性多步法的进一步扩展 在积分公式 （60）中，用个节点，的次插值多项式代替被积函数，得到的就是一般的线性多步法。如当时， （61）如当时， （62）特别地，用四个节点，的三次插值多项式代替被积函数，得到 （63）四个节点，的三次插值多项式代替，又得到 （64）其中误差分别为、，可见，公式（64）的精度比公式（63）的精度高。Adams 预报-校正公式

11、： （65）其中，说明：公式（65）自身不能提供最初所需的三个初值，需要Runge-Kutta 法的配合。4.4 一般线性多步方法一般的步线性多步方法可以写成 （8.24）这里，当时，上述公式就是显式公式，时为单步法.1、数值积分构造方法数值积分方法是构造线性多步法的一种途径.Adams方法就是利用插值多项式进行积分得出来的.又如，由Simpson求积公式有，其中，对应的数值公式. （8.25）称为Simpson方法.Simpson方法是四阶方法，它的局部截断误差为.2、Taylor展开构造方法构造线性多步方法的另一重要途径是利用Taylor展开方法.例如确定公式（8.24）中的个参数，方法是

12、将写成，并将所有项在进行Taylor展开，比较两边相同幂次前的系数，若， （8.26）则公式（8.24）是阶方法，局部截断误差主项为，其中.一般应用的线性多步方法都是大于或等于一阶的方法，这样的方法亦称为与微分方程（8.3）相容的方法. 阶方法的构造是选择参数，使（8.26）中的为.例6 推导最高阶的二步线性多步方法.解 1） 二步显式线性多步方法为.共四个待确定参数.由公式（8.26），令解得.又因为，则得阶数最高的二步显式线性多步方法为，其局部截断误差主项为，是三阶方法.2） 二步隐式线性多步方法为.共五个待确定参数.由公式（8.26），令解得.则得Simpson方法.【注】 两种多步法的构造公式中，数值积分方法是有局限的，它只对能将微分方程（8.3）转化为等价的积分方程的情形才能适用，而用Taylor展开则可构造任意多步法公式.3、几个常用的线性多步法的预测-校正公式（1） ilne预测校正法(ilnes Predic