-------------
第九章
常微分方程的数值解法
在自然科学的许多领域中，都会遇到常微分方程的求解问题。然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。这些方法可以给出解的近似表达式，通常称为近似解析方法。还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。利用计算机解微分方程主要使用数值方法。 我们考虑一阶常微分方程初值问题
?dy??f(x,y) ?dx??y(x0)?y0在区间[a, b]上的解，其中f (x, y)为x, y的已知函数，y0为给定的初始值，将上述问题的精确
解记为y(x)。数值方法的基本思想是：在解的存在区间上取n + 1个节点
a?x0?x1?x2???xn?b
这里差hi?xi?1?xi，i = 0,1, …, n称为由xi到xi+1的步长。这些hi可以不相等，但一般取
b?a。在这些节点上采用离散化方法，（通常用数值积分、微分。泰勒n展开等）将上述初值问题化成关于离散变量的相应问题。把这个相应问题的解yn作为y(xn)的近似值。这样求得的yn就是上述初值问题在节点xn上的数值解。一般说来，不同的离散化导致不同的方法。
§1 欧拉法与改进欧拉法
1.欧拉法 1．对常微分方程初始问题
成相等的，这时h?
?dy (9.1)??f(x,y) ?dx? (9.2)?y(x0)?y0
用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。 欧拉法是解初值问题的最简单的数值方法。从（9.2）式由于y (x0) = y0已给定，因而可以算出
y'(x0)?f(x0,y0)
设x1 = h充分小，则近似地有：
y(x1)?y(x0)?y'(x0)?f(x0,y0)
h (9.3)
-------------
-------------
记 yi?y(xi) i?0,1,?,n 从而我们可以取
y1?y0?hf(x0,y0)
作为y (x1)的近似值。利用y1及f (x1, y1)又可以算出y(x2)的近似值：
y2?y1?hf(x1,y1)
一般地，在任意点xn+1 = (n + 1)h处y(x)的近似值由下式给出
yn?1?yn?hf(xn,yn)
（9.4）
这就是欧拉法的计算公式，h称为步长。 不难看出，近似解的误差首先是由差商近似代替微商（见(9.3)）引起的，这种近似代替所产生的误差称为截断误差。还有一种误差称为舍入误差，这种误差是由于利用（9.4）进行计算时数值舍入引起的。 例9.1 用欧拉法求初值问题
0.9?y?y'??1?2x??x0?0?y(x0)?10.9y代入欧拉法计算公式。就得 1?2x0.9yn1?2xn???yn?n?0,1,?,5(9.5)(9.6)
当h = 0.02时在区间[0, 0.10]上的数值解。
解 把f(x,y)??yn?1?yn?h
?0.018??1??1?2xn?
具体计算结果如下表：
n xn yn y(xn) ?n = y(xn) - yn 0 0 1.0000 1.0000 0 1 0.02 0.9820 0.9825 0.0005 2 0.04 0.9650 0.9660 0.0005 3 0.06 0.9489 0.9503 0.0014 4 0.08 0.9336 0.9354 0.0018 5 0.10 0.9192 0.923 0.0021
在上表中y(xn)列，乃是初值问题（9.5）、（9.6）的真解
y(x)?(1?2x)?0.45
在xn上的值。?n为近似值yn的误差。从表中可以看出，随着n的增大，误差也在增大，所
-------------
-------------
以说，欧拉法计算简便，对一些问题有较大的使用价值，但是，它的误差较大，所得的数值解精确度不高。 2．改进欧拉法 为了构造比较精确的数值方法，我们从另一角度重新分析一下初值问题。一般说来，一阶方程的初值问题与积分方程
y(x)?y0??xx0f(t,y(t))dt （9.7）
是等价的，当x = x1时，
y(x1)?y0??x1x0f(t,y(t))dt （9.8）
要得到y(x1)的值，就必须计算出（9.8）式右端的积分。但积分式中含有未知函数，无法直接计算,只好借助于数值积分。假如用矩形法进行数值积分则
?x1x0f(t,y(t)dt?f(x0,y(x0))(x1?x0)
因此有
y(x1)?y0?f(x0,y(x0))(x1?x0)?y0?hf(x0,y0)
可见，用矩形法计算右端的积分与用欧拉法计出的结果完全相同。因此也可以说欧拉法的精度之所以很低是由于采用矩形法计算右端积分的结果。可以想象，用梯形公式计算（9.8）式右端的积分，可期望得到较高的精度。这时
1f(t,y(t))dt?{f(x0,y(x0))?f(x1,y(x1))}(x1?x0) x02将这个结果代入（9.3）并将其中的y(x1)用y1近似代替，则得
1 y1?y0?h[f(x0,y0)?f(x1,y1)]
2这里得到了一个含有y1的方程式，如果能从中解出y1，用它作为y(x1)的近似值，可以认为比用欧拉法得出的结果要好些。仿照求y1的方法，可以逐个地求出y2, y3，…。一般地当求出yn以后，要求yn+1，则可归结为解方程：
h yn?1?yn??f(xn,yn)?f(xn?1,yn?1)?
2这个方法称为梯形法则。用梯形法则求解，需要解含有yn+1的方程式，这常常很不容易。为此，在实际计算时，可将欧拉法与梯形法则相结合，计算公式为
?x1
(0)?yn?yn?hf(xn,yn)??1?(k?1)h(k)y?y?f(x,y)?f(x,y?n?1nnnn?1n?1)2???k?0,1,2,... （9.9）
(0)这就是先用欧拉法由(xn, yn)得出y(xn+1)的初始近似值yn然后用（9.9）中第二式进行迭代，?1，(k)(k?1)(k)反复改进这个近似值，直到yn为止，并把yn?1取作y(xn+1)?1?yn?1??（?为所允许的误差）
的近似值yn+1。这个方法就叫改进欧拉方法。
(0)(1)显然，应用改进欧拉法，如果序列yn?1,yn?1,?收敛，它的极限便满足方程
-------------