常微分方程的数值解法

第九章 常微分方程的数值解法在自然科学的许多领域中，都会遇到常微分方程的求解问题。

然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。

在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。

这些方法可以给出解的近似表达式，通常称为近似解析方法。

还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。

利用计算机解微分方程主要使用数值方法。

我们考虑一阶常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy在区间[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 可以不相等，但一般取成相等的，这时nab 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 hx 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 xy x f 219.0),(+-=代入欧拉法计算公式。

就得5,,1,021018.01219.01 =⎪⎪⎭⎫ ⎝⎛+-=+-=+n y xy x hy y n n nnn 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．改进欧拉法为了构造比较精确的数值方法，我们从另一角度重新分析一下初值问题。

一般说来，一阶方程的初值问题与积分方程⎰+=xx 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）式右端的积分，可期望得到较高的精度。

这时⎰-+≈10)())}(,())(,({21))(,(011100x x x x x y x f x y x f dt t y t f 将这个结果代入（9.3）并将其中的y (x 1)用y 1近似代替，则得)],(),([21110001y 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，则可归结为解方程：[]),(),(2111+++++=n n n n n n y x f y x f hy 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 收敛，它的极限便满足方程[]),(),(2111+++++=n n n n n n y x f y x f hy y 即序列的极限可取作y n +1。

可以证明，如果yf∂∂有界，则只要h 取得适当小，上述序列必定收敛。

这样当h 取得充分小，就可保证上述迭代过程收敛到一个解。

当步长h 取得适当时，欧拉方法算出的值已是较好的近似，因此改进欧拉法收敛很快，通常只需二、三次迭代即可。

如果迭代很多步仍不收敛，这表明表长h 选的过大，应缩小步长后再计算。

通常把（9.9）叫做预报校正公式，其中第一式叫预报公式，第二式叫校正公式。

这个公式还可写为⎪⎪⎪⎩⎪⎪⎪⎨⎧++==++=+),(),(2121121211k y h x hf k y x hf k k k y y n n n n n n (9.9)’3．公式的截断误差现在来考察两个公式的截断误差：11)(++-n n y x y 有多大？这里假定前一步得的结果y n =y (x n )是准确的。

写出y (x n +1)的泰勒展开式为+''++=+=+)(!2)(')()()(21n n n n n x y h x hy x y h x y x y 由欧拉法得)(')(),(1n n n n n n x hy x y y x hf y y +=+=+两式相减得)(0)(2)(2211h x y h y x y n n n =+''=-++即欧拉法的截断误差为0(h 2)，当h 0时它与h 2是同阶无穷小量。

对于改进的欧拉方法，我们以迭代一次的预报校正格式(9.9)’为例来说明。

因为)('))(,(),(1n n n n n x hy x y x hf y x hf k ===+''+=+⎭⎬⎫∂∂+⎩⎨⎧∂∂+=⎭⎬⎫⎩⎨⎧+∂∂+∂∂+=++=++=)()(')](,[)(')](,[)](,[)](,[)](,[)](,[])(,[),(221112n n n n n n n n n n n n n n n n n n n x y h x hy x y x f y x y x y x f x h x y x hf x y x f y k x y x f x h x y x f h k x y h x hf k y h x hf k 用它们代入（9.9）式第二式即得+''++=+''+++=+)(2)(')()(2)('2)('2221n n n n n n n n x y hx hy x y x y h x y h x y h y y这里末把含有h 的三次幂以上的项写出，因此有)(0)(311h y x y n n =-++即迭代一次的预报校正格式（9.9）的截断误差为0(h 3)。

可见改进的欧拉方法比欧拉法的阶提高了。

例9.2⎪⎩⎪⎨⎧=-=1)0(2y y x y dxdy 在区间[0, 1.5]上，取h = 0.1，求解。

解：（1）用欧拉法计算公式如下：1.0,1201==⎪⎪⎭⎫⎝⎛-+=+h y y x y h y y n n n n n（2）用迭代一次的改进欧拉法计算公式如下：⎪⎪⎭⎫ ⎝⎛-+=+n n n n y x y h y y 2)0(1[]1.0,1222),(),(20)0(11)0(1)0(111==⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛-+⎪⎪⎭⎫ ⎝⎛-+=++=++++++h y y x y y x y h y y x f y x f h y y n n n n n n n n n n n n n 本题的精确解为x x y 21)(+=，可用来检验近似解的精确程度。

计算结果如下表：x n欧拉法y n 迭代一次 改进欧拉法y n准确解n n x x y 21)(+=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.52.1200522.0142072.000000§2 龙格――库塔法由上节知道，截断误差的阶是衡量一个方法精度高低的主要依据。