常微分方程数值解法

1、第五章 常微分方程数值解法 1引言 在实际工作中常常会遇到求解常微分方程的定解问题。虽然我们已经学过一些常微分方程的定解问题的解法，但是那只是对一些特殊类型的方程给出了解析方法。而实际问题中归纳出来的常微分方程往往不能用解析方法来求解。有时虽然能求出其解析式子，但真要求它在某一点的值时，还会遇到一些麻烦。见下例。 例 已知初值问题 求其解。 易知，。若要求它在某 一点的值，还要求积分，而这个积分是不能求出的，需要用数值积分来求解。因此不如直接用数值解法求其解了。以下我们就研究常微分方程的数值解法。本章主要考虑下列一阶方程的定解问题（又称初值问题）。 （1-1）这里假设函数满足Lipschitz

2、（李卜西兹）条件。即满足：其中，L为某一常数。这是保证问题（1-1）有唯一解。下面就研究关于问题（1-1）的数值解法。所谓数值解法，就是直接寻求解函数在一系列离散点上的近似值，。其中，称为步长，今后如果不特殊说明，我们总是假定是等步长的。即有 ，此时节点，。由于为已知，所以自然设想利用这个已知信息求出的近似值，然后由求得的近似值，如此继续下去，这就是初值问题数值解法的一般思想。称为“步进法”。2 Euler方法（折线法）2-1 Euler公式这一方法是初值问题数值解中最简单的一个方法，其精度不高。所以实际计算中很少应用。但在某种程度上反映了数值方法的基本思想，且在此基础上得到的某些改进的方法目

3、前还在使用，因此我们有必要介绍一下这种方法。Euler公式的推导方法很多，在此我们用Tarlor展开的方法。 其中， ，由（1-1）式知， 当h充分小时，略去，即得 取近似，并写成等式得， （2-1）此即称为Euler公式。由于略去的是，可见这就是误差。但这只是在计算第n步时产生的误差，并非是从一开始到现在所产生的误差。因此我们称这个误差为局部截断误差。即在假定的假设下，计算时产生的误差，称为局部截断误差。由此知Euler公式的局部截断误差约为。2-2后退的Euler公式后退的Euler公式的推导同样也有许多方法。在此，我们用差商代替导数的方法。因为假设，又有（1-1）知，取近似整理得， （2

4、-2）此即称为后退的（隐式）Euler公式。公式（2-2）在计算时要用。这样的公式称为隐式公式，而象公式（2-1）的式子称为显式公式。对于隐式公式，在计算时要先给提供一个初值，然后再用给定的公式开始计算，通常称为迭代法。对（2-2）式的迭代公式如下， （2-3）直到为止。 以下研究式（2-2）的局部截断误差。 又 代入上式得，又有 下式减上式，且注意到已假设，故得， = 后退的Euler公式的局部截断误差约为。 2-3 梯形公式 今考察以上两个公式，可见它们的局部截断误差相差一个负号，即有 以上两式相加除以2，则可以消去，从而得到略去，则得一公式， （2-4）此称为梯形公式。其几何意义如图。

5、注意到（2-4）式为隐式公式，同样需要迭代法来求解。其迭代公式为， （2-5）以下分析迭代公式（2-5）的收敛性。（2-4）减去（2-5）得由于已假设满足Lipschitz（李卜西兹）条件。即满足：所以有，取，则有 从而收敛（带条件收敛）。 2-4 改进的 Euler公式 由（2-5）式知，用梯形公式计算时要反复求函数值，因而工作量是很大的。而改进的Euler公式是指在用梯形公式迭代时只迭代一次便取作，这就是， 预测： 校正： （2-6）以上称为预测-校正系统。也称为改进的 Euler公式为便于上机常采用以下形式， （2-7） 例 （见教材） 在以上推导后退的Euler公式时，我们用差商代替了

6、导数。但若改用中心差商代替导数就有 取近似并整理得，从而得， 用此公式作为预测值，用梯形公式作为校正值，便得到以下的预测-校正系统， 预测： 校正： （2-8）（2-8）与（2-6）相比，突出的特点是它们有相同的精度。不过（2-8）式是两步公式，而前面几式是单步公式。两步以上的公式称为多步法，多步法不能自开始（自起步），需要借助其他单步法才能开始进行。而单步法可以自开始。 3 Runge-Kutta方法 在1中，我们曾经用Taylor展开的方法推导了Euler公式。那时我们略去了项，自然想到如果多取几项，也就是略去项（p取的更大些）应该得到精度更高的公式。这就是Taylor级数法。它的一般公式

7、如下， （3-1）（3-1）式的局部截断误差为 = 一般地有以下定义， 定义 若一种方法的局部截断误差为，则称该方法具有p阶精度。 由（3-1）知，当p=1时，（3-1）式即为Euler公式。显然Euler公式具有1阶精度。同理，后退的Euler公式也具有1阶精度；而梯形公式具有2阶精度。由（3-1）知，p越大精度越高。但是它要用到f(x,y)的高阶导数，而求f(x,y)的高阶导数是不容易求出的。因此，Tarlor级数法在实际问题中很少应用，它的作用在于启发我们去探索更好的方法。3-1 Runge-Kutta方法的基本思想根据微分中值定理，存在，使得 由方程（1-1）知，有 （3-2）令 称为

8、区间上的平均斜率。由此可见，只要对这个平均斜率提供一种算法，由（3-2）式便可以得到一种计算公式。若取，即取点的斜率作为整个区间上的平均斜率，这时，则得Euler公式 若取，即取点的斜率作为整个区间上的平均斜率，这时，则得后退的Euler公式 若取，令 即取和两点斜率的算术平均作为整个区间上的平均斜率，则得下列公式， （3-3）这显然是改进的Euler公式（即梯形公式）。我们知道，改进的Euler公式是具有2阶精度，而Euler公式和后退的Euler公式只有1阶精度。可见用两个点斜率的算术平均作为平均斜率比只取一个点的斜率作为平均斜率精度要高些。由此得，如果设法在内多预测几个点的斜率值，然后将

9、它们加权平均作为平均斜率，则有可能构造出具有更高精度的计算公式。这就是Runge-Kutta方法的基本思想。3-2 二阶Runge-Kutta方法公式（3-3）是一个特殊的二阶Runge-Kutta方法。它用了，这两个点上的斜率值。我们现在在内取 ， 将这两个点上的斜率值，的线性组合作为平均斜率。显然，而，对于其中的我们用Euler公式来计算得， （3-4）为此有， （3-5）其中含有三个待定参数。确定这三个待定参数的原则是使得公式（3-5）具有二阶精度。为此将在点作二元函数Taylor展开，得 将，代入（3-5）的第一式得 而二阶Taylor级数法公式为 比较二者系数得到，三个待定参数应满足

10、 （3-6）满足（3-6）形如（3-5）的所有公式均称为二阶Runge-Kutta公式。特别当，p=1时就是改进的Euler公式。3-3三阶Runge-Kutta方法为提高精度，我们现在在内取 ， 三个点，用这三个点上的斜率值，的线性组合作为平均斜率。显然，而的预测值我们用来计算。从而得， （3-7）为了确定（3-7）中的待定参数。我们同样用将，作Taylor展开，然后代入（3-7）中的。通过与三阶Taylor级数法比较系数，则得待定参数满足的条件为, （3-8）以及 （3-9）参数满足（3-8）（3-9）的形如（3-7）的公式统称为三阶Runge-Kutta公式。它也是一族公式。用同样的方法

11、我们还可以得到四阶Runge-Kutta公式。Runge-Kutta公式的一般表示为 取r=4,用同样的方法即可以得到四阶Runge-Kutta公式。下面就是一个经典的四阶Runge-Kutta公式。 （3-10）可以证明其截断误差为。应注意的是， Runge-Kutta公式的推导基于Taylor展开方法，所以它要求解函数有很好的光滑性。即y(x)要具有所要求的导数。若不然，用高阶Runge-Kutta公式可能不如用低阶Runge-Kutta公式效果好。4 单步法的收敛性与稳定性4-1 收敛性定义2 若一种数值方法对于任意固定的，单步法 （4-1）满足 ，则称该方法是收敛的。例 考虑 （4-2

12、）解 （4-2）的Euler公式为，所以有 ，因为， 。对于一般单步法 （4-1）有以下结论，定理 假设单步法有p阶精度，且增量函数关于y满足Lipschitz条件。即对任意，总存在常数L0,使 （4-3）则单步法（4-1）收敛。且整体截断误差为。定理的证明从略。作为应用我们考虑梯形公式（改进的Euler方法）的收敛性。梯形公式的增量函数为 而以上两式相减，并注意到不等式的性质，得 又由于我们已假设函数f(x,y)关于y满足Lipschitz条件，所以上式变为 整理得， 所以，梯形公式的增量函数关于y满足Lipschitz条件。从而收敛。4-2 稳定性定义3 若一种数值方法在节点上的值有大小的

13、扰动，而以后节点(mn)上的值所产生的扰动不超过，则称该方法是稳定的。由于稳定性的讨论比较复杂，所以通常对模型方程 （） （4-4）进行讨论，以得到稳定性条件。所以稳定性的定义还可以如下描述，定义4 用单步法（4-1）解模型方程（4-4）所得到的稳定性方程，若，则称此法是绝对稳定的。由，所得区间称为稳定区间。为在节点值上的扰动值。今考虑用Euler公式解模型方程（4-4）时的稳定区间。为此设为理论值，为实际计算值，则有 理论值 实际计算值 两式相减得稳定性方程， 其中为第n次扰动。显然，要保证稳定的条件是，即有。解之得Euler公式的稳定区间为对于后退的Euler公式，则有两式相减得稳定性方程

14、，解得由于，对任意步长h0总有，所以我们通常说后退的Euler公式是无条件稳定的。或者说后退的Euler公式的稳定区间为。同样方法我们可以讨论其他公式的稳定区间。比如，改进的Euler公式（2-6）。由于公式为预测： 校正： 对于模型（4-4）变为 预测： 校正： 即 两式相减得稳定性方程， 得稳定性条件为 解之得稳定区间为其他公式不再讨论。 例 见教材 5 线性多步法以上研究的都是单步法，这一节我们研究线性多步法。一般线性多步法的公式为 （5-1）其中，。当，（5-1）为显式公式；当，（5-1）为隐式公式。适当选择（5-1）式中的待定参数，就可以得到一系列的线性多步法公式。确定待定参数时常用

15、的方法有数值积分法和Taylor展开法，在此我们用Taylor展开法来确定待定参数，因为这种方法更具有一般性，且更容易理解。 设 分别在点作Taylor展开，得， 以上两式代入（5-1）并整理得 要使上式具有p阶精度，即局部截断误差为。又由的Taylor展开式 知，需要符合到项，所以应有下式成立 （5-2）当（5-2）式成立时，局部截断误差为 （5-3） 以下我们具体的考虑两个四步方法。先考虑显式公式 要使上式有四阶精度，根据（5-2）式知其系数应满足 （5-4）因为共有7个待定参数，而只有5个方程，所以（5-4）的解是不唯一的。若取，则得 。此时的公式为 （5-5）其中，。（5-5）称为四步

16、Adams显式公式。再考虑隐式公式（用代替），即 要使上式有四阶精度，根据（5-2）式知其系数应满足 （5-6）同样有7个待定参数，而只有5个方程，所以（5-6）的解是不唯一的。若取，则得。此时的公式为 （5-7）其中，。（5-7）称为四步Adams隐式公式。用以上的方法还可以得到Milne公式以及Hamming公式等。作为以上方法的应用我们考虑下例，例 解初值问题用显式二步公式。 （5-8）试确定待定参数，使方法阶数尽可能高；并求局部截断误差。解 由于以上公式不易记忆，这里我们直接用Taylor展开方法。因为 其中，代入（5-8）式得，考虑Taylor级数法公式为使方法阶数尽量高，需有解之得，。此时公式为三阶。所得二步法为 将，代入得 而Taylor级数法公式当p=3时相应的项为局部截断误差为 作为练习请你用以上方法解下例： 例 解初值问题用显式二步公式。 试确定待定参数，使方法阶数尽可能高；并求局部截断误差的主项。 例 证明线性多步法 存在的一个值，使方法是四阶的。解 将展