第九章 常微分方程的数值解法

第九章 常微分方程的数值解法

在自然科学的许多领域中，都会遇到常微分方程的求解问题。然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。这些方法可以给出解的近似表达式，通常称为近似解析方法。还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。利用计算机解微分方程主要使用数值方法。
我们考虑一阶常微分方程初值问题

在区间[a, b]上的解，其中f (x, y)为x, y的已知函数，y_0为给定的初始值，将上述问题的精确解记为_y(x)。数值方法的基本思想是：在解的存在区间上取n + 1个节点

这里差，i = 0,1, …, n称为由x到x+1的步长。这些h可以不相等，但一般取成相等的，这时。在这些节点上采用离散化方法，（通常用数值积分、微分。泰勒展开等）将上述初值问题化成关于离散变量的相应问题。把这个相应问题的解y作为y(x)的近似值。这样求得的y就是上述初值问题在节点x上的数值解。一般说来，不同的离散化导致不同的方法。

§1 欧拉法与改进欧拉法

1.欧拉法
1．对常微分方程初始问题


用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。从（9.2）式由于y (x) = y已给定，因而可以算出

设x_1 = _h充分小，则近似地有：
(9.3)
记
从而我们可以取

作为y (x)的近似值。利用y_1及_f (x, y)又可以算出y(x)的近似值：

一般地，在任意点x+1 = (n + 1)h处y(x)的近似值由下式给出
（9.4）
这就是欧拉法的计算公式，h称为步长。
不难看出，近似解的误差首先是由差商近似代替微商（见(9.3)）引起的，这种近似代替所产生的误差称为截断误差。还有一种误差称为舍入误差，这种误差是由于利用（9.4）进行计算时数值舍入引起的。
例9.1 用欧拉法求初值问题

当h = 0.02时在区间[0, 0.10]上的数值解。
解 把代入欧拉法计算公式。就得

具体计算结果如下表：

在上表中y(x)列，乃是初值问题（9.5）、（9.6）的真解

在x上的值。为近似值y的误差。从表中可以看出，随着n的增大，误差也在增大，所以说，欧拉法计算简便，对一些问题有较大的使用价值，但是，它的误差较大，所得的数值解精确度不高。
2．改进欧拉法
为了构造比较精确的数值方法，我们从另一角度重新分析一下初值问题。一般说来，一阶方程的初值问题与积分方程
（9.7）
是等价的，当x = x_1时，
（9.8）
要得到_y(x)的值，就必须计算出（9.8）式右端的积分。但积分式中含有未知函数，无法直接计算,只好借助于数值积分。假如用矩形法进行数值积分则

因此有

可见，用矩形法计算右端的积分与用欧拉法计出的结果完全相同。因此也可以说欧拉法的精度之所以很低是由于采用矩形法计算右端积分的结果。可以想象，用梯形公式计算（9.8）式右端的积分，可期望得到较高的精度。这时

将这个结果代入（9.3）并将其中的y(x)用y_1近似代替，则得

这里得到了一个含有_y_1的方程式，如果能从中解出_y_1，用它作为_y(x)的近似值，可以认为比用欧拉法得出的结果要好些。仿照求y_1的方法，可以逐个地求出y, _y，…。一般地当求出y以后，要求y+1，则可归结为解方程：

这个方法称为梯形法则。用梯形法则求解，需要解含有y+1的方程式，这常常很不容易。为此，在实际计算时，可将欧拉法与梯形法则相结合，计算公式为
（9.9）
这就是先用欧拉法由(x, y)得出y(x)的初始近似值，然后用（9.9）中第二式进行迭代，反复改进这个近似值，直到（e为所允许的误差）为止，并把取作y(x)的近似值y+1。这个方法就叫改进欧拉方法。
显然，应用改进欧拉法，如果序列收敛，它的极限便满足方程

即序列的极限可取作y+1。可以证明，如果有界，则只要h取得适当小，上述序列必定收敛。这样当h取得充分小，就可保证上述迭代过程收敛到一个解。当步长h取得适当时，欧拉方法算出的值已是较好的近似，因此改进欧拉法收敛很快，通常只需二、三次迭代即可。如果迭代很多步仍不收敛，这表明表长h选的过大，应缩小步长后再计算。通常把（9.9）叫做预报校正公式，其中第一式叫预报公式，第二式叫校正公式。这个公式还可写为
(9.9)’
3．公式的截断误差
现在来考察两个公式的截断误差：有多大？这里假定前一步得的结果y=y(x)是准确的。
写出y(x)的泰勒展开式为

由欧拉法得

两式相减得

即欧拉法的截断误差为0(h)，当h ®0时它与h_2是同阶无穷小量。
对于改进的欧拉方法，我们以迭代一次的预报校正格式(9.9)’为例来说明。因为


用它们代入（9.9）式第二式即得

这里末把含有_h的三次幂以上的项写出，因此有

即迭代一次的预报校正格式（9.9）的截断误差为0(h)。可见改进的欧拉方法比欧拉法的阶提高了。
例9.2

在区间[0, 1.5]上，取h = 0.1，求解。
解：
（1）用欧拉法计算公式如下：

（2）用迭代一次的改进欧拉法计算公式如下：


本题的精确解为，可用来检验近似解的精确程度。计算结果如下表：

x	欧拉法y	迭代一次 改进欧拉法y	准确解
0	1	1	1
0.1	1.1	1.095909	1.095445
0.2	1.191818	1.184096	1.183216
0.3	1.277438	1.260201	1.264911
0.4	1.358213	1.343360	1.341641
0.5	1.435133	1.416102	1.414214
0.6	1.508966	1.482956	1.483240
0.7	1.580338	1.552515	1.549193
0.8	1.649783	1.616476	1.612452
0.9	1.717779	1.678168	1.673320
1.0	1.784770	1.737869	1.732051
1.1	1.85118	1.795822	1.788854
1.2	1.917464	1.852242	1.843909
1.3	1.984046	1.907323	1.897367
1.4	2.051404	1.961253	1.949359
1.5	2.120052	2.014207	2.000000

§2 龙格――库塔法

由上节知道，截断误差的阶是衡量一个方法精度高低的主要依据。能否用提高截断误差阶来提高方法的精确度呢？回答是肯定的。本节介绍的泰勒级数法和龙格――库塔法就是基于这种思想构造出来的。
1．泰勒级数法
如果初值问题
（9.10）
的精确解y(x)在[x, x]上存在k + 1阶导数且连续，那么由泰勒公式
（9.11）
其中截断误差为
（9.12）
略去截断误差，并用近似值代替真值则得
（9.13）
用公式(9.13)解初值问题的数值方法称为泰勒级数法，当 k = 3时，（9.13）变为
（9.14）
这时的截断误差是
（9.15）
从截断误差的表示式中看出，如果微分方程（9.10）的真解y = y(x)为次数不超过3的多项式时，公式（9.14）精确地成立，因此（9.14）是3阶方法。
例9.3 导出用三阶泰勒级数法解方程

的计算公式
解:
因




故

而

其中表示f(x, y)对x的k阶偏导数在x = x点上的值。
泰勒级数法只要初值问题的真解充分光滑，就可以获得精确度较高的数值解。但是须计算y(x)的各阶导数，这当f (x, y)的表达式复杂时是很繁琐的。因此泰勒级数法一般只用于求“表头”（即开头几点的数值解，如y_1, _y, y, , y等）。另外，用上述级数法计算表头时，还可以得到选择步长h的信息。假定我们要求计算误差不超过e ，那么，当h满足条件
（A）
（B）
时，应该认为是最好的。因为，当条件（A）不满足时，达不到指定精确度，而当条件（B）不满足则表明h过小。
能否构造一种格式，既保留泰勒级数法精确度较高的优点，又避免过多的计算f (x, y)的各阶偏导数呢？下面介绍的龙格――库塔方法就能办到这一点。
2．龙格――库塔法
从理论上讲，只要函数y = y(x)在区间[a, b]上充分光滑，那么它的各阶导数值y(k)(x)与函数y(x)在区间[a, b]上某些点的值就相互有联系，就是说，函数值可用各阶导数值近似地表示出来，反之，各阶导数值也可用函数在一些点上值的线性组合近似地表示出来。事实上，欧拉法和改进欧拉法也可以看成是导数值用函数值的线性组合表示的特例，例如，改进欧拉法可以写成
（9.16）
此公式也可称为二阶龙格――库塔公式。
为了导出龙格――库塔法的一般公式，我们取如下的线性组合形式，
（9.17）
其中
（9.18）
即

而w_1, _w,…, w；b_1 =0, _b, b,…, b；a_21, _a, …, a除b_1=0外均为待定系数。适当选取这些系数，使得局部截断误差的阶尽可能高即可。
显然，当_g = _1时，（9.12）式就是欧拉公式。
下面我们先导出_g =_2时的公式。将_k_1, _k在同一点（x, y）泰勒展开，则有

（9.19）
将(9.19)代入(9.17)并与y(x+h)在x点的泰勒展开式：

逐项比较，令h、h_2项的系数相等，便得到

把_b_2作为自由参数来确定_w_1和_w_2，如取_b_2 = 1，则_w_1 = _w = ，a_21 = 1，这时（9.17）正好就是改进的欧拉方法，截断误差的阶为0(_h)。
对于g = _3的情形，我们也可以完全仿上述方法推导出三阶龙格――库塔公式。这时参数满足下列条件
（9.20）

（9.20）比较简单的一组解为：
_b =，b_3 = 1，_a_21 = ，_a_31 = -1，_a_32 = 2，_w_1 = ，_w_2­ = ，_w_3 =
将它代入（9.17）得
(9.21)
这就是三阶龙格――库塔公式。当然，参数的不同选取公式（9.21）就有不一样的形式，但它们的截断误差阶都是0(_h)。
通常人们所说的龙格――库塔法是指四阶而言的。我们可以仿照二阶的情形推导出此公式，不过太繁杂，此处从略，常用的四阶公式是
（9.22）
公式(9.22)的截断误差阶为0(h)。
龙格――库塔法有精确度高、收敛、稳定（在一定的条件下）计算过程中可以改变步长等优点，但仍需计算f (x, y)在一些点的值，如四阶龙格――库塔法每计算一步需要算四次f (x, y)的值，这就给实际计算带来一定的复杂性。因此，它与泰勒级数法一样，一般用于计算“表头”。
例9.4 用龙格――库塔法解初值问题
y’ = x_2 – y (0≤_x≤1) y(0) = 1 （9.23）
解 ： 取 h = 0.1，由（9.22）得

把初始条件x_0 = 0，_y_0 = 1，代入,得_k_1 = -1，_k_2 = -0.9475，_k_3 = -0.9501，_k_4 = 0.8950，将这些_k值代（9.22），得

重复上述步骤可算出y_2，_y_3，…，_y_10等。

§3 线性多步法

前面介绍的方法，统称为单步法，就是在计算_y+1时，只用到前面一步y的值。本节我们将介绍利用前边已经算出来的若干个值y，y+1，…，y-1，y，来求得y+1的高精度公式――线性多步方法。
我们已经知道初值问题
y’ = f (x, y) y (x) = y（9.24）
与积分方程
（9.25）
等价。本节介绍的方法，其基本思想是用一个插值多项式P(x)来代替（9.25）中的被积函数f (x, y)，然后用
（9.26）
代替(9.25)。这种方法实际上要分两步来做：
（1）求出开头几个点上的近似值，即计算“表头”；
（2）利用(9.26)逐步求后面点x上的值y。
选取不同的插值点作插值多项式，就会得出不同的数值解法，我们这里只讨论其中的一种，叫阿当姆斯（Adams）方法。
（一）阿当姆斯外推公式
为了简单起见，我们以k = 2为例，导出阿当姆斯外推公式。于是以x-2，x-1，x为节点作牛顿向后插值多项式P_2(_x)。
（9.27）
其中

而余项为
（9.28）
将（9.27）代入（9.26）式并经变量替换x = x + th，便得
（9.29）
在上述公式中被插值点x≤x≤x+1不包含在插值节点所决定的最大区间（x-2, x）内，故（9.29）称阿当姆斯外推公式(或称外插公式)，显然，它是显式的且每前进一步只计算一次f (x, y)的值即可。
公式（9.26）的截断误差，可由（9.28）及积分第二中值定理得

上式表明阿当姆斯外推法(当k = 2时)的截断误差的阶为0(h)。
同理，对k = 3的情形即可求得外推公式
（9.31）
而余项为
（9.32）
公式（9.29）和（9.31）要输入一个差分表。为此，我们将差分表示成函数值的和的形式：
（9.33）
其中是组合数，即

这样（9.29）可改写为
（9.34）

而（9.31）可改写为
（9.35）
公式（9.34）和（9.35）是常用的外推公式。
例9.5 求y’ = x + y, y (0) = 1
当x = 0.1到0.5，步长h = 0.1时的数值解。
解 先用前面讲过的方法计算出表头
y = 1； y = 1.11034；
y_2 = 1.24281； _y = 1.39972
将上述值代入（9.35）式，计算得：
y = 1.58364； y = 1.79742
（二）阿当姆斯内插法
根据插值理论知道，节点的选择对于精度有直接的影响。同样次数的内插公式比外插公式更精确。
如果（9.26）中的被积函数是以x-1，x，x+1为插值节点的内插多项式，即
（9.36）
将（9.36）代入（9.26）便得
（9.37）
把它改写成便于在计算机上实现的形式，就有
（9.38）
上式便是k = 1时的阿当姆斯内插公式，它的截断误差的阶为
（9.39）
同理，对k = 2的情形可导出公式
（9.40）
或
（9.41）
而 （9.42）
这就是常用的阿当姆斯内插公式。
比较外推法和内插法的截断误差可知，在同样利用k + 1个已知值时，阿当姆斯内插法的截断误差阶比外推法高一阶，因而内插法更精确。
但内插法（9.38）及（9.41）是隐式的，需要解方程，通常用迭代法求解。如用外推法（9.34）算出的值作为初始近似，然后相应地用内插法公式（9.41）进行迭代，即

（9.43）
若将（9.41）记成的形式，容易算出

因此，当h充分小时，可有

即迭代程序（9.43）收敛，并且h越小，收敛越快。当h选得适当时，一般迭代二、三次即可得到满足精度要求的y+1值。
同理可导出一般的阿当姆斯内插公式。
为了提高精度，经常把阿当姆斯外推法与内插联合起来交替使用。例如


第一个方程的精度为0(h)，用第二个方程迭代一次精度仍达0(h)。这样两个方程交替使用可达较好的效果。第一个方程是预报方程，第二个方程是较正方程。
例9．6 对例9.5用内插法求解
解 先用前面讲过的方法计算出表头
y = 1； y = 1.11034； y = 1.24281
再用（9.43）的第一个式子计算出，最后用（9.43）的第二个式子进行迭代，得
y = 1.39972
同样，可算出y_4及_y_5。

§4 解二阶常微分方程边值问题的差分法

考虑常微分方程的边值问题：
（9.44）
其中_p(x)，q(x)和f (x)均为[a, b]上给定的函数，a，b为已知数。
边值问题的主要特点是在两个端点上给出了定解条件。在以下讨论中我们总假定p(x)、q(x)及f (x)均为[a, b]上充分光滑的函数，且q(x)≤0，这时，边值问题（9.44）存在连续可微的解，且唯一。
用差分法解边值问题的主要步骤是：
（1）将区间[a, b]离散化；
（2）在这些节点上，将导数差商化，从而把微分方程化为差分方程；
（3）解差分方程――实际上就是解线代数方程组。
今将[a, b]区间用节点

分成N等分，其中x_0 = _a与x = b称为边界点，而x_1, _x, …, x称为内点。
再将y²(x)，y’(x)在内点x_1处的值，_y’(x)可表示为中心差商
（9.45）



从而y² 的中心差商为
（9.46）
略去上式中截断误差的阶为0(h)，并以近似值y代替真值y(x)，便得


同一、二阶中心差商分别代替（9.44）中的一、二阶导数，则得逼近（9.44）的差分方程
（9.47）
y = __a，y = b，i = 1，…，N – _1
其中_p= p(x)，q = q(x)，f = f(x)。经过简单整理，（9.47）可改写成
（9.48）
差分方程（9.48）实际上是含有N – _1个未知数，_y_1,…, _y的线代数方程组。
最后，解差分方程（9.48）就可得数值解y。为此，令
，
（9.49）
则有
（9.50）
这个方程组的系数矩阵是三对角矩阵，可以用追赶法求解。
例9.7 试用差分法解方程
（9.51）
解 将[0, 1]划分为四等分，即取，得五个节点

（9.51）的差分方程为
（9.52）
将它改写成
（9.53）
在每个内点列方程得
（9.54）
由追赶法递推公式解得：
y = 1.4855 y = 1.2802 y = 0.7753

习题
1．取步长h = 0.1用改进欧拉法解初值问题

试将计算结果与准确解相比较。
2．取步长h = 0.2再用四阶龙格――库塔方法解初值

并用前题比较结果。
3．取步长h = 0.2用四阶龙格――库塔方法解

4．下列各题先用龙格――库塔法求表头，然后用阿当姆斯法继续求以后各值
（1）
（2）
5．用差分法求方程

的数值解（h = 0.2）
6．求方程

的数值解（取h = 0.2）。