常微分方程的数值解法

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

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

⎪⎩⎪⎨⎧==0

0)()

,(y x y y x f dx

dy

在区间[a, b]上的解，其中f (x, y )为x, y 的已知函数，y 0为给定的初始值，将上述问题

的精确解记为y (x )。数值方法的基本思想是：在解的存在区间上取n + 1个节点

b x x x x a n =<<<<= 210

这里差i i i x x h -=+1，i = 0,1, …, n 称为由x i 到x i +1的步长。这些h i 可以不相等，但一

般取成相等的，这时n

a

b h -=。在这些节点上采用离散化方法，（通常用数值积分、微分。泰勒展开等）将上述初值问题化成关于离散变量的相应问题。把这个相应问题的解y n 作为y (x n )的近似值。这样求得的y n 就是上述初值问题在节点x n 上的数值解。一般说来，不同的离散化导致不同的方法。

§1 欧拉法与改进欧拉法

1.欧拉法

1．对常微分方程初始问题

(9.2)

)((9.1)

),(00⎪⎩⎪⎨

⎧==y x y y x f dx dy

用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。

欧拉法是解初值问题的最简单的数值方法。从（9.2）式由于y (x 0) = y 0已给定，因而可以算出

),()('000y x f x y =

设x 1 = h 充分小，则近似地有：

),()(')

()(00001y x f x y h

x y x y =≈-

(9.3)

记 ,n ,,i x y y i i 10 )(== 从而我们可以取

),(0001y x hf y y ==

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

),(1112y x hf y y +=

一般地，在任意点x n +1 = (n + 1)h 处y (x )的近似值由下式给出

),(1n n n n y x hf y y +=+

（9.4）

这就是欧拉法的计算公式，h 称为步长。

不难看出，近似解的误差首先是由差商近似代替微商（见(9.3)）引起的，这种近似代替所产生的误差称为截断误差。还有一种误差称为舍入误差，这种误差是由于利用（9.4）进行计算时数值舍入引起的。 例9.1 用欧拉法求初值问题

⎪⎩⎪⎨⎧

==+-=)

6.9(01

)()5.9(219.0'00x x y y x y

当h = 0.02时在区间[0, 0.10]上的数值解。

解 把y x

y x f 219

.0),(+-

=代入欧拉法计算公式。就得

5

,,1,021018.01219

.01 =⎪⎪⎭

⎫ ⎝

⎛+-=+-=+n y x

y x h

y y n n n

n

n n

具体计算结果如下表：

n x n y n y (x n ) n = y (x n ) - y n 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 (x n )列，乃是初值问题（9.5）、（9.6）的真解

45.0)21()(-+=x x y

在x n 上的值。n ε为近似值y n 的误差。从表中可以看出，随着n 的增大，误差也在增大，所

以说，欧拉法计算简便，对一些问题有较大的使用价值，但是，它的误差较大，所得的数值解精确度不高。

2．改进欧拉法

为了构造比较精确的数值方法，我们从另一角度重新分析一下初值问题。一般说来，一阶方程的初值问题与积分方程

⎰

+

=x

x dt t y t f y x y 0

))(,()(0

（9.7）

是等价的，当x = x 1时，

⎰

+

=1

))(,()(01x x dt t y t f y x y （9.8）

要得到y (x 1)的值，就必须计算出（9.8）式右端的积分。但积分式中含有未知函数，无法直

接计算,只好借助于数值积分。假如用矩形法进行数值积分则

⎰

-≈1

)())(,()(,(0100x x x x x y x f dt t y t f

因此有

)

,()

))((,()(000010001y x hf y x x x y x f y x y +=-+≈

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

⎰

-+≈1

0)())}(,())(,({2

1

))(,(011100x x x x x y x f x y x f dt t y t f 将这个结果代入（9.3）并将其中的y (x 1)用y 1近似代替，则得

)],(),([2

1

110001y x f y x f h y y ++=

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

[]),(),(2

111+++++=n n n n n n y x f y x f h

y y

这个方法称为梯形法则。用梯形法则求解，需要解含有y n +1的方程式，这常常很不容易。为此，在实际计算时，可将欧拉法与梯形法则相结合，计算公式为

[]

⎪⎩

⎪⎨⎧=++=+=+++++,...

2,1,0)

,(),(2),()(11)1(1)0(1k y x f y x f h y y y x hf y y k n n n n n k n n n n n （9.9）

这就是先用欧拉法由(x n , y n )得出y (x n +1)的初始近似值)

0(1+n y ，然后用（9.9）中第二式进行

迭代，反复改进这个近似值，直到ε<-+++)(1)1(1k n k n y y （为所允许的误差）为止，并把)

(1k n y +取

作y (x n +1)的近似值y n +1。这个方法就叫改进欧拉方法。

显然，应用改进欧拉法，如果序列 ,,)

1(1)0(1++n n y y 收敛，它的极限便满足方程