计算方法第7微分方程数值解

计算方法北京工业大学应用数理学院杨中华第七章常微分方程数值解

有很多实际问题是用微分方程建立其数学模型，但是目前只能对少量典型的微分方程求出相应的解析形式的解，而对于绝大多数的微分方程问题，很难或者根本不可能得到它的解析解。因此，研究求解微分方程的数值解是非常有意义的。

我们知道微分方程是关于未知函数及导数的方程，包括常微分方程和偏微分方程。常微分方程所涉及的未知函数是一元函数；如果未知函数是多元的，则称为偏微分方程，偏微分方程数值解有另外独立的课程介绍。本章仅讨论常微分方程初值问题的数值解方法。7.1

常微分方程的初值问题

考虑如下常微分方程初值问题

所谓求初值问题的数值解就是求其解析解

在点

处的近似值

一般采用等距节点的方式，也就是取

，称为步长，并令

，从

求得

，再从

求得

，依次求解就可得到y=y(x)在

各点的数值解

。

(7.1)(7.2)

定义7.1如果确定数值解

的方法仅使用

则称该方法为求解初值问题的单步法；如果确定

的方法需要

则称为该方法求解初值问题的多步法。

一般的单步法形如

在(7.3-1)式中右端无

项可直接计算，称为显式单步法；而(7.3-2)公式无法直接计算

必须通过解方程得到，称为隐式单步法。初值问题(7.3)(7.3-1)(7.3-2)如

求解微分方程的线性多步法形如

例如

这两个计算公式均是线性二步法，而(7.4-1)式可以直接计算

，是显式线性二步法；注意到(7.4-2)式中右端因此(7.4-2)是一个关于

的方程，是隐式线性多步法。

初值问题(7.4)(7.4-1)(7.4-2)

定义7.2设y(x)是初值问题的解析解，

是按某数值方法在

点所求的计算解，称

为该数值方法的整体截断误差或总体截断误差；称

为该数值方法在

点的局部截断误差。

显然局部截断误差是从

到

一个算法步骤所产生的误差，也就是在

准确的前提下，

的误差；初值问题(7.5)(7.6)

所谓总体截断误差是算法从

计算

、再从

计算

、…、直到从

计算

这

n+1个步骤所积累的总误差。

因此局部截断误差与总体截断误差是不同的但有关联，特别如果在计算

之前的

是准确时，局部截断误差

与总体截断误差

是相同的。

求解微分方程初值问题数值方法的收敛性问题是假定每一步计算都没有舍入误差的前提下，讨论步长

时，方法的总体截断误差是否趋于零的问题。初值问题

定义7.3

如果数值方法的局部截断误差可表示为

则称这种数值方法的阶数是p。

数值方法的阶数决定了该方法的收敛性和收敛速度。

定义7.4

设

是初值问题在

点的两个近似初值，

是某数值方法从这两个初值开始所求的近似解，如果存在常数

，当步长满足

且

时，有下式成立

则称该数值方法是稳定的。初值问题

定义7.4其含义是如果求解微分方程的数值方法以

的任意两个初值作为起点，所得两近似解的误差将被两个初值的差所界定，进而可以推出，稳定的数值方法在求解过程中具有自动校正的性质，也就是从

开始求解的近似解序列

会落入积分曲线

一个邻近通道内。如下图所示：初值问题

定理7.1设

在区域

内连续且关于y满足Lipschitz条件：存在常数L，使

对任意

和任意

成立，则初值问题(7.1)存在连续可微且唯一的解。

此处列出定理7.1的原因是因为以后的很多结论均是以Lipschitz条件为前提的。初值问题(7.7)7.2Euler方法

1.Euler公式

先从研究一阶微分方程初值问题的几何直观入手，问题(7.1)式可以表示为平面问题：

这表明该微分方程的解析解

所表示的积分曲线上每一点

切线斜率等于函数

的值，也就是说，积分曲线上每一点的切线方向可用二元函数

表示，如图。

求解此微分方程初值问题的困难在于二元函数

中的y是所求的未知函数，除初始点的

外其他点函数值均未知，对于其他点的函数近似值，可使用如下方法求解：

初值条件

作为积分曲线的出发点

以曲线在

处的切线为方向，前进到直线

得到

，即

欧拉方法

这里

称为步长。同样以

点处的

值为方向，前进到直线

到

，即

如图所示，再从

得到

，从

得到

，…

按此过程我们得到如下迭代公式：

这就是著名的Euler公式。 1)Euler公式是在

处用向前差商推导而得； 2)Euler公式的本质是用当前点的斜率来预测曲线下一个点

的位置，然后用折线来近似所求的积分曲线； 3)当h很小时折线对积分曲线的偏离不会太大，而且Euler 公式具有自动校正的特点，即折线会向积分曲线靠拢欧拉方法(7.8)

例7.1

用Euler公式求解初值问题

在

处的近似值，取步长

。

解：

由Euler公式

详细计算结果列表如下

该问题的解析解为

，因此

在

处，

，3位有效数字n10.021.0200000020.041.0408000030.061.0624160040.081.0848643250.101.10816161欧拉方法

2.截断误差

现在我们对Euler方法进行误差分析，也就是讨论Euler公式的截断误差。1)局部截断误差

设

是微分方程初值问题的准确解，令

则

为Euler公式在

处的局部截断误差，考虑

在

处的Taylor展开式

将

带入此式并注意到

是

的解析解，可得欧拉方法

从而局部截断误差为

其中

称为误差的主项，从(7.9)式可以看出局部截断误差与

同价。2)总体截断误差

设

是由Euler公式得到的，称

为Euler公式在

处的总体截断误差，记作

若

关于y满足Lipschitiz条件，即存在正数L，使得对一切

有

则Euler公式的总体截断误差为欧拉方法(7.9)(7.10)

证明：设

是Euler公式在

处由准确的y(x)计算得到的结果，则有

其中第一部分是局部截断误差有

考虑第二部分，由Lipschitiz条件

将两部分综合考虑则有欧拉方法

注意到

，当

时有

再注意到

则

有界，即总体截断误差与h同阶，也就说明

也就说明Euler方法是收敛的，且为一阶的方法。欧拉方法

3.后退Euler公式

前述讨论的Euler公式是对

用向前差商近似导数得到的，如果使用向后差商，则有

由此导出另一种形式的Euler公式，称为后退Euler公式

后退Euler公式与(向前)Euler公式的局部截断误差是同阶的，都是

，但后退Euler公式与Euler公式有着本质的区别。Euler公式是关于

的一个直接计算公式，这类公式称为显式；而后退Euler公式的右端含有未知的

，它实际上是关于

的一个方程，这类公式为隐式。欧拉方法(7.11)

通常选用第二章的迭代法求解该隐式方程，而只需迭代一、二步就可以达到精度要求。

向前Euler公式仅仅根据

一个点的信息来构造

，而后退Euler公式是根据

两个点的信息来构造

，可以预期后退Euler公式将优于向前Euler公式，实际计算表明，后退Euler公式更稳定。

4.梯形公式和改进Euler公式1)梯形公式

前述的Euler公式是用直观方式导出的，其实Euler公式也可以用积分方式导出，加以改进可以得到另一个更有效的求解微分方程的公式，以下是导出过程。欧拉方法

设

是初值问题的解，因此有

，由此得到

在此公式推导过程中，是用左矩形的面积近似定积分的计算，显然定积分的计算如果使用梯形积分公式会更精确些，由此得到

称

为求解常微分方程的梯形公式。欧拉方法(7.11)

2)梯形公式的局部截断误差

设

是初值问题的准确解，记

为梯形公式在

处的局部截断误差，考虑

在

处的Taylor展开式

从第二个Taylor展开式可得到

将其代入第一个Taylor展开式可得到欧拉方法

其局部截断误差为

从而局部截断误差与h3同阶，其误差主项为

仿照Euler方法总体截断误差的推导过程，可以证明梯形法的总体截断误差为

，即梯形法是一个二阶方法。欧拉方法(7.13)

3)改进Euler公式

梯形公式的局部截断误差要比Euler公式高一阶，提高了精度，但因为梯形公式是隐式的，必须通过求解方程得方式得到

，算法相对复杂计算量也大。为了简化算法，通常在使用求解非线性方程的迭代法时只迭代一两次就转入下一步计算，简化算法如下：

先用Euler公式求得一个初步的近似值

，称为预测值，预测值的精度可能较差，再用梯形公式将其校正一次得

称为校正值，迭代格式如下：

该迭代格式称为预测—校正公式，也称为改进Euler公式。欧拉方法(7.14)

例7.2用改进Euler公式求解例7.1中的初值问题。

解：

预测

校正

详细计算结果列表如下

n10.021.020000001.0204000020.041.041208001.0416160830.061.063248401.0636647240.081.086138021.0865627550.101.109894011.11032732欧拉方法

该问题的解析解为

，因此

在

处，

，有5位有效数字，比Euler公式计算结果多2位有效数字。欧拉方法

4)改进Euler公式的局部截断误差

虽然改进Euler公式是简化的梯形公式，仅用一次迭代作为所求隐式方程的解，但是所得算法的局部截断误差却与

同阶，因此改进Euler法也是二阶的方法。证明如下：

设y=y(x)是初值问题的准确解，若记

考虑改进Euler公式一步产生的误差，即局部截断误差，

令

是梯形公式的局部截断误差，

是改进Euler公式的局部截断误差，则欧拉方法

即改进Euler公式局部截断误差与

同阶，从而可以知道改进Euler法是二阶的方法。欧拉方法7.3Runge-Kutta法

1.Runge—Kutta法的基本思想

前节给出的Euler公式和改进Euler公式的局部截断误差分别与

和

同阶，显然两公式所得计算解的精度后者将优于前者，本节先探讨两者优劣的原因，再引导我们构造出更为有效的计算公式。

考虑在

上对y=y(x)应用微分中值定理

令

，上式可写成如下形式

则Euler公式可表示为(7.15)(7.16)

如果把

看作是

相对于

的校正量(通常是未知的)，则Euler方法实际上是用一个点

处的信息(包括步长和积分曲线在此点上的斜率)构造

作为校正量

的近似值。

改进Euler公式可以表示为

而改进Euler方法是用两个点

、

处的信息构造

和

的算术平均值作为

的近似值，也就是说改进Euler方法用了更多的信息构造校正量

的近似值，当然如此改进的Euler方法效果应该优于Euler方法。龙格|库塔方法(7.17)

再仔细观察，在改进Euler公式中，把

点处的值

作为预测值代入

计算出

，然后再用两值的平均值去近似

。进一步，如果能在区间

内多预测几个点的

值，并用它们的加权平均来作为

的近似值，即

可以预期会构造出具有更高精度的计算公式，这就是Runge—Kutta法的基本思想。龙格|库塔方法(7.18)

2．二阶Runge—Kutta法

考虑用

和

的加权平均值来近似

的一般形式，设

要使其局部截断误差达到

，则其参数必须满足以下方程组

注意到方程组有三个方程，4个未知量，因此方程组有无穷多个解，不同的解将对应不同的计算公式，通常有如下三种参数的取法。龙格|库塔方法(7.19) 1)改进Euler法

取

，就是改进Euler公式。

2)Heun公式

取

，则得到Heun公式。

3)修正的Euler公式或中点公式

取

，所得公式称为修正的Euler公式或中点公式，即龙格|库塔方法(7.20)(7.21)

我们可以通过中点公式的几何意义来体会二阶Euler公式为何更为有效。首先考察积分曲线在

点附近是凸函数的情况，如图：龙格|库塔方法

从图上看出Euler公式的计算值

对应C点，真正的积分曲线在

点的准确值

对应A点，中点公式的计算值

对应B点附近。

事实上，当积分曲线是凸函数且h足够小时，将有龙格|库塔方法

其中

就是图中三角形与

重合的一条直角边的边长，该边长较之于

更接近于准确的

，从而经此校正之后的计算解

对应的B点更接近于准确解的A点。

对于积分曲线在

点附近是凹函数的情况，也容易看出中点公式优于Euler公式，如下图所示。

3.四阶Runge—Kutta法

按照Runge—Kutta法的基本思想，构造四阶Runge—Kutta法也就是用

的加权平均值来近似

，因此令

欲使局部截断误差达到

，即

，类似于二阶Runge—Kutta法，经过一系列较为复杂的推导过程，得到有13个参数的11个方程组(见文献[7])。由于方程的个数少于未知量的个数，因此方程组有无穷多个解，此处仅给出两个常用的四阶公式。龙格|库塔方法(7.22) (1)标准四阶Runge—Kutta公式龙格|库塔方法(7.23) (2)Gill公式(7.24)

例7.3用标准四阶Runge—Kutta公式求解例7.1中的初值问题

解：

与准确解比较高达7位有效数字。龙格|库塔方法

4．变步长Runge—Kutta法

在实际计算中，步长h的选择是一个十分重要的问题。若步长h取得太小，将增加许多不必要的计算量；若步长

取得太大,其计算精度很难得到保证。因此，仿照数值积分方法，常微分方程的数值解也有一个选择步长问题。可使用如下方法：

设

是从

点以h为步长使用四阶Runge—Kutta公式计算得到的，

是从

点以

为步长两次使用四阶Runge—Kutta法计算过程得到的。对于给定精度

，如果相邻两次计算满足

则停止计算，以最后计算结果为

的近似值；否则步长折半继续计算直到满足精度要求为止。龙格|库塔方法(7.25)7.4线性多步法

所谓多步法就是在逐步迭代的求解过程中，计算

时已经算出近似值

，如果能够充分利用前面多步的信息来预测

，则可能获得较高的精度，这就是所谓的线性多步法的基本思想。

1.线性多步法的一般公式

利用前面得到的多步信息

来计算

最常用的方法是线性多步法，其一般公式为

其中

为常数，

为

的近似值，

，当

时该式为显式，否则为隐式。多步法的本质就是用

的线性组合来构造

。(7.26)

定义7.2

设y(x)是微分方程初值问题的准确解，如果多步法公式的局部截断误差满足

则称该多步法具有p阶精度，或称方法是p阶的。显然p愈大结果愈好。 2.线性二步法

设有如下3个线性二步法公式

前两个公式是具有2阶精度的显式线性二步法，第三个公式是利用Simpson求积公式求得具有4阶精度的隐式线性二步法。线性多步法(7.27)

3.Adams外推公式 Adams外推公式利用

的信息来构造计算

的公式，其形式如下

我们首先给出

的(7.28)式(当

时是Euler公式)，注意到

，将

在

点Taylor展开

将其代入到(7.28)式，则有

如果我们将

在

点Taylor展开线性多步法(7.28)(7.29)(7.30)

欲使(7.28)式的精度足够高，应使

中的p尽可能的大，也就是使(7.29)和(7.30)两式相同的项数尽可能的多，显然应取

则局部截断误差

，由此得到具有二阶精度的Adams公式：线性多步法(7.31)

Adams外推公式是一种显式格式,所以也称为Adams显式公式，仿照上述方法我们还可以推导出

情况的Adams公式，此处不再赁述(最常用的是

情况4阶精度的Adams公式)。线性多步法(7.