数值分析第9章 常微分方程数值解

1、许多实际问题的数学模型是微分方程或微分方程的定解问题。如物体运动、电路振荡、化学反映及生物群体的变化等。常微分方程可分为线性、非线性、高阶方程与方程组等类；线性方程包含于非线性类中，高阶方程可化为一阶方程组。若方程组中的所有未知量视作一个向量，则方程组可写成向量形式的单个方程。因此研究一阶微分方程的初值问题 (9-1)的数值解法具有典型性。常微分方程的解能用初等函数、特殊函数或它们的级数与积分表达的很少。用解析方法只能求出线性常系数等特殊类型的方程的解。对非线性方程来说，解析方法一般是无能为力的，即使某些解具有解析表达式，这个表达式也可能非常复杂而不便计算。因此研究微分方程的数值解法是非常必要

2、的。只有保证问题(9-1)的解存在唯一的前提下，研究其数值解法或者说寻求其数值解才有意义。由常微分方程的理论知，如果(9-1)中的满足条件（1）在区域上连续；（2）在上关于满足Lipschitz条件，即存在常数，使得则初值问题(9-1)在区间上存在惟一的连续解。在下面的讨论中，我们总假定方程满足以上两个条件。所谓数值解法，就是求问题(9-1)的解在若干点处的近似值的方法。称为问题(9-1)的数值解，称为由到的步长。今后如无特别说明，我们总假定步长为常量。建立数值解法，首先要将微分方程离散化，一般采用以下几种方法：(1) 用差商近似导数在问题(9-1)中，若用向前差商代替，则得用其近似值代替，所

3、得结果作为的近似值，记为，则有这样，问题(9-1)的近似解可通过求解下述问题(9-2)得到，按式(9-2)由初值经过步迭代，可逐次算出。此方程称为差分方程。需要说明的是，用不同的差商近似导数，将得到不同的计算公式。(2) 用数值积分法将问题(9-1)中的微分方程在区间上两边积分，可得(9-3)用,分别代替,，若对右端积分采用取左端点的矩形公式，即同样可得出显式公式(9-2)。类似地，对右端积分采用其它数值积分方法，又可得到不同的计算公式。(3) 用Taylor多项式近似。把在点处Taylor展开，取一次多项式近似，则得设，略去余项，并以代替，便得以上三种方法都是将微分方程离散化的常用方法，每一

4、类方法又可导出不同形式的计算公式。其中Taylor展开法，不仅可以得到求数值解的公式，而且容易估计截断误差。上面我们给出了求解初值问题(9-1)的一种最简单的数值公式(9-2)。虽然它的精度比较低，实践中很少采用，但它的导出过程能较清楚地说明构造数值解公式的基本思想，且几何意义明确，因此它在理论上仍占有一定的地位。1简单的数值方法和基本概念1.1 Euler法与向后Euler法一、Euler法Euler方法就是用差分方程初值问题(9-4)的解来近似微分方程初值问题(9-1)的解，即由公式(9-4)依次算出的近似值。从几何上看，微分方程在平面上确定了一个向量场：点处的方向斜率为。问题(9-1)的

5、解代表一条过点的曲线，称为积分曲线，且此曲线上每点的切向都与向量场在这点的方向一致。从点出发，以为斜率作一直线段，与直线交于点，显然有，再从出发，以为斜率作直线段推进到上一点，其余类推，这样得到解曲线的一条近似曲线，它就是折线。因此Euler方法又称为Euler折线法。二、向后Euler法在微分方程离散化时，用向后差商代替导数，即，则得到如下差分方程(9-5)用这组公式求问题(9-1)的数值解称为向后Euler法。向后Euler法与Euler法形式上相似，但实际计算时却复杂得多。Euler法计算的公式中不含有，这样的公式称为显式公式；向后Euler法计算的公式中含有，称为隐式公式。显式公式与隐

6、式公式各有特点。显式公式的优点是使用方便，计算简单，效率高。其缺点是计算精度低，稳定性差；隐式公式正好与它相反，它具有计算精度高，稳定性好等优点，但求解过程很复杂，一般采用迭代法。为了结合各自的优点，通常将显式公式与隐式公式配合使用，由显式公式提供迭代初值，再经隐式公式迭代校正。上面隐式公式中，在求解时，为已知，是方程的根。一般说来，这是一个非线性方程，因此我们通过构造简单迭代法来求解。迭代格式为由于满足Lipschitz条件，所以由此可知，只要，迭代法就收敛到解。1.2 梯形公式利用数值积分方法将微分方程离散化时，若用梯形公式计算式(9-3)中右端积分，即并用代替，则得计算公式(9-6)这就

7、是求解初值问题(9-1)的梯形公式。梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为由于函数关于满足Lipschitz条件，所以其中为Lipschitz常数。因此，当时，迭代法就收敛到解。1.3 局部截断误差与方法的精度为了刻画近似解的准确程度，引入局部截断误差与方法精度的概念。定义 假设在某一步的近似解是准确的，即（这个假设称为局部化假设）。在此前提下，用某公式推算所得，我们称为该公式（即该方法）的局部截断误差。定义9.2 如果某种方法的局部截断误差是则称该方法是阶方法，或具有阶精度。显然越大，方法的精度越高。1）Euler法的截断误差假设问题的解充分光滑，且前步计算结果是精确的，即，于

8、是Euler法的截断误差是(9-7)这里称为局部截断误差主项。显然。2）向后Euler法的截断误差。计算公式是将对用微分中值定理，有（在与之间）将在处Taylor展开于是将方程的解作Taylor展开因此故(9-8)3）梯形法的计算公式是将对用微分中值定理，有（在与之间）将在处Taylor展开于是将方程的解作Taylor展开因此故(9-9)所以梯形法是二阶方法。1.4 改进的Euler法我们看到，虽然梯形方法提高了精度，但其算法复杂，在应用迭代公式进行实际计算时，每迭代一次，都要重新计算函数的值，而迭代又要反复进行若干次，计算量很大。为了控制计算量，通常只迭代一两次就转入下一步的计算，这就简化了

9、算法。具体地说，我们先用欧拉公式求得一个初步的近似值，称之为预测值，预测值的精度可能很差，再用梯形公式将它校正一次得，称为校正值。这样的预测校正系统通常称为改进的欧拉公式。即预测：校正：为了编制程序上机，将上式改写成(9-10)算法(1) 输入，整数，初值(2) 置 输出(3) 计算 输出(4) 若，置，转3；否则停机。改进欧拉法的截断误差。计算公式是在处作Taylor展开式，注意到，有于是(9-11)所以改进Euler法是二阶方法。2 龙格-库塔（Runge-Kutta）法2.1 Runge-Kutta法的基本思想及一般形式设初值问题(9-1)的解，按微分中值定理，必存在，使(9-12)图9

10、-3其中叫做在上的平均斜率。只要对平均斜率提供一种算法，公式(9-12)就给出一种数值解公式。例如，用代替，就得到前进Euler公式；用替代，就得到后退Euler公式；如果用的算术平均值替代，则可得到2阶精度的梯形公式，如图9-3。可以设想，如果在上能多预测几个点的斜率值，用它们的加权平均值替代，就有望得到具有较高精度的数值解公式，这就是Runge-Kutta法的基本思想。Runge-Kutta公式的一般形式是(9-13)其中是在点的斜率预测值。均为常数。选取这些常数的原则是使公式(9-13)具有尽可能高的精度。公式(9-13)叫做级显式Runge-Kutta公式，简称RK方法。2.2 二阶R

11、K法的推导的RK方法的近似公式为(9-14)这里均为待定常数，我们希望适当选取这些系数，使公式阶数尽量高。先展开，按照二元函数Taylor级数：，得(9-15)为了叙述方便，把及其偏导数中的省略不写，将式(9-15)代入(9-14)的第一式得再展开，注意到公式得于是，局部截断误差为(9-16)要使式(9-16)的局部截断误差为，则应要求的系数为零，于是有(9-17)方程有4个未知数，3个方程，所以有无穷多组解，它的每组解代入式(9-14)得到的近似公式，局部截断误差均为，故这些方法统称为二阶方法。例如，取，得，近似公式为(9-18)这就是改进Euler公式。如果取，有，近似公式为(9-19)这

12、也是常用的二阶公式，称为中点公式。注：对于一般函数，由于，所以不论参数如何选取，只能有。这说明(9-14)至多是二阶的方法。类似地，对和的情形，通过更复杂的计算，可以导出三阶和四阶RK公式，其中常用的三阶和四阶RK公式为(9-20)和(9-21)式(9-21)称为经典的四阶RK公式，通常说四阶RK方法就是指用式(9-21)求解。算法(1) 输入(2) 置(3) 计算输出(4) 若，则置转3；否则停机。例4 用Euler法（），改进Euler法（）和经典R-K方法（）计算初值问题.其结果如下表9-1表9-1Euler方法()改进Euler法()经典R-K法()00000注：1）通过比较RK四阶经

13、典公式和前面用Euler法、改进Euler法（它是一种二阶RK方法）的计算结果，我们发现，四阶RK方法的精度高得多。就计算量来说，虽然Euler法、改进Euler法每步只需计算一个或二个函数值，而四阶RK方法每步需计算四个函数值，但由于放大了步长，计算量几乎相同。可见，选择有效的算法是非常重要的。2）RK方法的导出基于Taylor展开，故它要求问题的解具有较高的光滑度。当解充分光滑时，四阶RK方法确实优于改进Euler法。如果解的光滑性较差，则用四阶RK方法求得的数值解，其精度可能反而比改进Euler法的差。因此，在实际计算时，需根据问题的具体情况选择合适的算法。3）当时，RK公式的最高阶数恰

14、是，当时，RK公式的最高阶数不是，如时，RK公式的最高阶数仍为4，时，RK公式的最高阶数仍为5。由此看来，四阶RK公式是最经济的。2.3 变步长的RK方法若问题(9-1)的解函数变化是不均匀的，用等步长求数值解，则可能产生有些点处精度过高，有些点处精度过低的情况，为保证一定的精度，必须取较小的步长。这样做既增加了计算量，也导致误差的严重积累，因而实际计算时，往往采用事后估计误差、自动调整步长的RK方法。设用阶方法计算，记从出发，以步长计算一步所得的近似值，以步长计算两步得的近似值。由于阶公式的局部截断误差为，且在小区间上变化不大，故有(9-22)(9-23)将式(9-24)乘以减(9-23)得

15、从而有(9-24)这样，可以事后误差估计(9-24)中的大小来选择合适的步长。变步长的RK方法与数值积分做法相同，可以从的值来选择合适的步长，从而得到合乎精度要求的。具体作法是：设要求精度是，如果，就反复加倍步长进行计算，直到为止。这时上一次的步长就是合适的步长，上一次计算所得的结果，就是合乎精度要求的；如果，则反复折半步长，直到为止。这时，最后一次步长就是合适的步长，而这最后一次的计算结果就是满足精度要求的。这种计算过程中自动选择步长的方法，叫变步长方法。以上所介绍的各种数值解法都是单步法，这时因为它们在计算时，都只用到前一步的值，单步法的一般形式是 (9-25)其中为增量函数，例如Eule

16、r方法的增量函数为，改进Euler法的增量函数为3 相容性、收敛性与稳定性3.1 相容性与收敛性上面所介绍的方法都是用离散化的手段，将微分方程初值问题化为差分方程初值问题求解的。这些转化是否合理？即当时，差分方程是否能无限逼近微分方程，差分方程的解是否能无限逼近微分方程初值问题的准确解，这就是相容性与收敛性问题。用单步法(9-25)求解初值问题(9-1)，即用差分方程初值问题 (9-26)的解作为问题(9-1)的近似解，如果近似是合理的，则应有 (9-27)（体会：要使得差分方程能无限逼近微分方程，我们将微分方程的精确解代人差分方程中，当时，应该满足）其中为问题(9-1)的精确解。因为故由(9

17、-27)得如果增量函数关于连续，则有 (9-28)如果单步法的增量函数满足条件(9-28)，则称单步法(9-25)与初值问题(9-1)相容。通常称(9-28)为单步法的相容条件。满足相容条件(9-28)是可以用单步法求解初值问题(9-1)的必要条件。容易验证Euler法和改进Euler法均满足相容性条件。事实上，对Euler法，增量函数为自然满足条件(9-28)，改进Euler法的增量函数为因为连续，从而有所以Euler法和改进Euler法均与初值问题(9-1)相容。一般地，如果单步法有阶精度（），则其截断误差为上式两端同除以，得令，如果连续，则有所以的单步法均与问题(9-1)相容。由此即得各

18、阶RK方法与初值问题(9-1)相容。定义9.4 一种数值方法称为是收敛的，如果对于任意初值及任意固定的，都有其中为初值问题(9-1)的精确解。如果我们取消局部化假定，使用某单步法公式，从出发，一步一步地推算到处的近似值。若不计各步的舍人误差，而每步都有局部截断误差，这些局部截断的积累就是整体截断误差。定义9.5我们称为某数值方法的整体误差。其中为初值问题(9-1)的精确解，为不计舍人误差时用某数值方法从开始，逐步得到的在处的近似值（不考虑舍人误差的情况下，截断误差的积累）。定理设单步法(9-25)具有阶精度，其增量函数关于满足Lipschitz条件，问题(9-1)的初值是精确的，即，则单步法的

19、整体截断误差为证明 由已知，关于满足Lipschitz条件，故存在，使得对任意的及，都有记，因为单步法具有阶精度，故存在，使得从而有反复递推得因为，即，又，于是所以推论1 设单步法具有()阶精度，增量函数在区域：上连续，且关于满足Lipschitz条件，则单步法是收敛的。当在区域上连续，且关于满足Lipschitz条件时，改进Euler方法，各阶RK方法的增量函数在区域上连续，且关于满足Lipschitz条件，因而它们都是收敛的。关于单步法收敛的一般结果是：定理设增量函数在区域上连续，且关于满足Lipschitz条件，则单步法收敛的充分必要条件是相容性条件(9-28)。3.2 稳定性稳定性与收

20、敛性是两个不同的概念，收敛性是在假定每一步计算都准确的前提下，讨论当步长时，方法的整体截断误差是否趋于零的问题。而稳定性则是讨论舍人误差的积累能否对计算结果有严重影响的问题。 若一种数值方法在节点值上有一个大小为的扰动，于以后各节点上产生的偏差均不超过，则称该方法是稳定的。我们以Euler方法为例进行讨论。假设由于舍人误差，实际得到的不是而是，其中是误差。由此再计算一步，得到把它与不考虑舍人误差的Euler计算公式相减，并记，就有其中。如果满足条件， (9-29)则从到的计算，误差是不增的，可以认为计算是稳定的。如果条件(9-29)不满足，则每步误差将增大。当时，显然条件(9-29)不可能满足

21、，我们认为问题本身具有先天的不稳定性。当时，为了满足收敛性要求(9-29)，有时要很小，这样步数就相当多，这时误差的累积可能是十分严重的，出现数值不稳定现象。现在我们考虑的初值问题 (9-30)其解是。今设在的初值有误差，实际求解的问题是 (9-31)它的准确解是。因此如果用很精确的方法，求出(9-31)的解，则它和(9-30)的解在处误差为，而在处的误差则是。它随着的增大而增大，与所选的数值方法无关，是问题本身固有的特性。再看一个的例子，初值问题为其准确解是。用欧拉法求解有 (9-32)若取，则欧拉公式的具体形式为同样讨论有初始误差，则在的误差是。数值不稳定。式(9-32)中取，则欧拉公式的

22、具体形式为。对初始误差，在的误差是。数值稳定。以上讨论表明稳定性不但与方法有关，也与步长的大小有关，当然也与方程中的有关。为了只考察数值方法本身，一般只检验数值方法用于求解模型方程的稳定性，模型方程为 (9-33)其中为复数。这个方程的分析比较简单，一般的方程可以通过局部线性化为这种形式，例如在的邻域略去高阶项，再作变量替换就得到的形式。事实上，方程可简写为，作变换即可得。对于个方程的方程组，可以线性化为，其中为的Jacobi矩阵。若有个不同的特征值，则可对角化为，再作变换，得到个非耦合的方程，其中可以复数，所以一般讨论(9-33)中的为复数。对于方程(9-33)，若，类似以上分析，认为方程是

23、不稳定的。所以我们考虑的情形，这时不同的数值方法可能是数值稳定或不稳定的。当一个单步法用于试验方程，从计算一步得 (9-34)其中依赖于所选的方法。因为通过点试验方程的解曲线（它满足）为，而一个阶单步法的局部截断误差在时有，所以有 (9-35)这样可以看出是的一个近似值。由(9-34)可以看到，若计算中有误差，则计算时将产生误差，所以有下面定义。定义9.6 如果(9-34)式中，则称单步法(9-25)是绝对稳定的。在复平面上复变量满足的区域，称为方法(9-25)的绝对稳定区域，它与实轴的交称为绝对稳定区间。在上述定义中，规定严格不等式成立，是为了和线性多步法的绝对稳定性定义一致。事实上，时也可

24、以认为误差不增长。(1) Euler方法的稳定性Euler方法用于模型方程(9-33)，得，所以有。所以绝对稳定条件是，它的绝对稳定区域是复平面上以为中心的单位圆。而为实数时，绝对稳定区间是。(2) 梯形公式的稳定性对模型方程，梯形公式的具体表达式，即，所以梯形公式的绝对稳定条件为。化简得，因此梯形公式的绝对稳定区域为平面的左平面。特别地，当为负实数时，对任意的，梯形公式都是稳定的。图9-4 图9-5(3) RK方法的稳定性与前面的讨论相仿，将RK方法用于模型方程(9-33)，可得二、三、四阶RK方法的绝对稳定区域分别为当为实数时，三、四阶显式RK方法的绝对稳定区域分别为，。表9-2 常用几个

25、公式的绝对稳定区间方法阶数区间Euler公式（显式）1后退Euler公式（隐式）1梯形公式（隐式）2二级R-K公式（显式）2二级R-K公式（隐式）2三级R-K公式（显式）3经典R-K公式（显式）4例5 分别取和0.2，用经典R-K方法计算解 本例，分别为和，前者在绝对稳定区间内，后者则不在，这里列出计算误差表9-3115251256253125从表9-3可以看出，当时，虽然计算结果的相对误差很大，但计算过程是稳定的。而当时，计算过程不稳定。这是因为时，属于稳定区间，当时，不属于稳定区间。对于一般方程，在考虑稳定性对步长的限制时，常用代替（因为前面变换中），只要步长的选取使得在所用方法的绝对稳定

26、区域内即可。例如若用Euler方法求解初值问题因为且，要使，只要4 线性多步法前面所讲的各种数值解法，都是单步法。因为它们在计算时，都仅仅用到前面的一个信息；如果在计算值时，能够比较充分地利用前面的已知信息，如，那么就可望使所得到的更加精确。这就是多步法的基本思想。多步法中最常用的是线性多步法，其一般公式是 (9-36)其中均为常数，。当时，上式称为显式，否则称为隐式。线性多步法与Runge-Kutta法都是高精度方法，不同的是线性多步法是利用前面已求得的在各节点处的近似值及导数近似值的线性组合来逼近的，而Runge-Kutta法则是用在内若干点处的一阶导数预测值的线性组合，来逼近在区间内的平

27、均斜率，从而求得的近似值。构造线性多步法的公式有多种途径，最常用的是Taylor展开方法和数值积分的方法.定义设是初值问题(9-1)的精确解，线性多步法(9-36)在上的局部截断误差为若，则称方法(9-36)是阶的，则称方法(9-36)与问题(9-1)的第一个方程相容的。4.1 线性多步法公式的推导利用Taylor展开导出线性多步公式(9-36)的基本方法是：将线性多步公式(9-36)在处进行Taylor展开，然后与在处的Taylor展开式相比较，要求它们前面的项重合，由此确定参数。设初值问题(9-1)的解充分光滑，记，把它在处展成Taylor级数则有 (9-37)及 (9-38)用替代式(9

28、-37)、(9-38)中的，则得及把这两个式子代入式(9-36)，得(9-39)为了使式(9-36)具有阶精度，只须使式(9-39)的前项与在处的Taylor展开式(9-40)的前项系数对应相等即可。对比关于的同次项系数，得到确定的方程组： (9-41)可见，只要式(9-36)的系数满足式(9-41)，方法(9-36)就具有阶精度。式(9-40)减去(9-39)，可得局部截断误差为 (9-42)4.2常用的线性多步法公式式(9-41)共有个方程，个待定常数：，只要，即，就可以由式(9-41)，定出具有阶精度的公式(9-36)的系数。下面介绍几个常用的线性多步法公式。一、阿达姆斯（Adams）公

29、式取，有(9-43)此时方程有9个未知数，5个方程，有四个自由变量。特别取，可解得相应的线性多步法公式为(9-44)因为，式(9-44)称为Adams显式公式，它是四阶公式，局部截断误差为 (9-45)如果令，由式(9-43)可得解相应的线性多步法公式为(9-46)因为，式(9-46)称为Adams隐式公式，它是四阶公式，局部截断误差为 (9-47)二、米尔恩（Milne）公式对方程组(9-42)，令，解出相应的线性多步法公式 (9-48)称为Milne公式，其局部截断误差为 (9-49)Milne公式是四阶四步显式公式。三、海明（Hamming）公式对方程组(9-41)，取，并令，可得Ham

30、ming公式 (9-50)其局部截断误差为 (9-51)Hamming公式是四阶三步隐式公式。由于线性多步法公式(9-36)，只有在知道了前面的之后，才能使用。所以开头的的值必须先用同阶的单步法（Taylor级数法或Runge-Kutta法）求出，然后才能用线性多步法公式（例如，若使用四阶线性多步法公式，必须用同阶单步法，算出）。例6 分别用四阶Adams显式和隐式公式求初值问题的数值解，取。解 根据题意，由四阶Adams显式公式(9-43)有由四阶Adams隐式公式(9-44)有由上式可解出利用精确解求出起步值（一般可用四阶RK公式计算起步值）后，按上面的公式计算，结果见下表表9.4Adam

31、s显式法Adams隐式法00从表9.4可以看出，四阶Adams隐式法比显式法的精度高，比较这两种方法的局部截断误差公式(9-45)与(9-47)，可以说明这一现象。一般地，同阶的隐式法比显式法精确，而且数值稳定性也好。但在隐式公式中，通常很难用迭代法求解，这样又增加了计算量。因此实际计算时，很少单独使用显式公式或隐式公式，而是将它们联合使用：先用显式公式求出的预测值，记作，再用隐式公式对预测值进行校正，求出的近似值。4.3 预测-校正系统一般地，采用同阶的显式公式与隐式公式配对使用，即把由显式求出的（记）作为的预测值，然后再代入隐式公式进行校正，求出更接近的值。这样就构成了预测-校正系统。采用的预测校正系统有两种：Adams显式-隐式公式 (9-52)Milne-Hamming预测校正系统 (9-53)说明：(1) 以上两种预测校正公式均为四阶公式，其起步值通常用四阶RK公式计算。(2) 有时为提高精度，校正公式可迭代进行多次，但迭代次数一般超过3次。为减少一次迭代所产生的误差，常常用局部截断误差进一步修正预测值与校正值