常微分方程数值解法的误差分析 (2)

2、常微分方程初值问题描述

在自然科学和经济的许多领域中，常常会遇到一阶常微分方程的初值问题

?dy??f?x,y?,a?x?b (1) ?dx?y?x0??y0?这里f?x,y? 是充分光滑，即关于x或y满足李普希茨条件的二元函数，y0是给定的初值，y?x0??y0称为初始条件。

3、数值解法的基本思想与途径

一阶微分方程的初值问题(1)的解y?x?是区间?a,b?上的连续变量x的函数，因而问题(1)实际上是一个连续性的问题，求这个问题的数值解，就是要求在区间

?a,b?上的若干个离散点处的函数近似值，例如：

a?x0?x1?...?xn?b，

然后计算出解y?x?的近似值

y?x0?,y?x1?,...,y?xn?.

一般常取x0,x1,...,xn为等距离的点，即

x1?x0?x2?x1?...?xn?xn?1?h

?0,1,.. .n,或 xi?a?i,hi称h为步长。

建立数值方法的第1步，就是把连续性问题(1)通过一定的方法化为在给定的

n?1个点上的近似的差分方程的初值问题，称这个过程为离散化。常用离散化的方法如下：

(一)用差商替代导数

在点xi处的导数yi??xi?可以近似地表示成差商 y?yiy??xi??i?1,

h从而把初值问题(1)化为差分问题

- 2 -

?yi?1?yi?f?xi,yi?,?hi?0,1,...n??y?x0??y0,? (2)

其中yi表示解y?x?在点xi处的近似解，即yi?y?xi?。

当然，用差商来近似地表示导数，方法不是唯一的，这里所用的是所谓的向前差商。

(二)Taylor展开法

在一点(例如点xi)的附近，y?x?的同次数的近似多项式中的Taylor多项式

hp?p?y?xi?h??y?xi??hy??xi??...?y?xi?p!

为最好。其中p为一正整数。通过微分方程y??f?x,y?，便可以逐次把各阶导数y?,y??,?在xi处的值表示出来。

(三)数值积分法

对微分方程y??f?x,y?在区间?xi,xi?1?上求积分，得

y?xi?1??y?xi???xi?1xif?x,y?x??dx,?, i?0,1于是，初值问题(1)便可以近似地化为

?y?y?xi?1f?x,y?x??dx,i?i?1?xii?0,1,?n. ???y?x0??y0,这样，关于上式右端的积分，可以用数值积分方法计算其近似值。

4、数值解的分类

常微分方程初值问题的数值解法一般分为两大类：

单步法：所谓单步法是指这类方法在计算yn?1时，只用到前一步的值

xn?1,xn,yn,然后逐步往下计算。这个算法的代表是龙格---库塔算法，简称R—K

方法。四阶显示Runge-Kutta方法是求解普通常微分方程初值问题数值解法中的重要方法，而隐式Runge-Kutta公式是求解刚性常微分方程初值问题的重要方法。

多步法：这类方法在计算yn?1时，除了用到前一步的值xn?1,xn,yn,之外，还要用到

xn?p，yn?p?p?1,2,?,k;k?0?

- 3 -

这前面k步的值，这个算法的代表就是阿达姆斯(Adams)方法。

5、问题(1)解的存在惟一性定理

一个常微分方程是不是有特解呢?如果有，又有几个呢?这是微分方程论中一个基本的问题，数学家把它归纳成基本定理，叫做存在和唯一性定理。因为如果没有解，而我们要去求解，那是没有意义的；如果有解而又不是唯一的，那又不好确定。因此，存在和唯一性定理对于微分方程的求解是十分重要的。这个重要的存在和唯一性就是下面列出的著名的存在惟一性定理。

???y????中连定理?1? 如果f?x,y?在带形区域R???x，y?a?x?b，续，且关于y满足Lipchiz条件：即存在正常数L，使得

f续可微解

对所有的x??a,b?以及任何

?x,1y???f,x2?y?L1?y2

,yy?y?x?。

y1，y2都成立，那么初值问题(1)存在惟一的连

二、几种常用的数值解法及其误差分析

首先来了解下初值问题数值解中最简单的一种方法——欧拉法。

1、单步法

Euler折线法发生的历史背景。在微分方程研究之初，瑞士数学家L．Euler(1707.4-1783.9)做出了开创性的工作。他和其他一些数学家在解决力学、物理学问题的过程中创立了微分方程这门学科。在常微分方程方面，Euler在1743年发表的论文中，用代换

y?ekx扩给出了任意阶常系数线性微分方程的古典解

法，最早引入了“通解”和“特解’’的概念。

1768年，Euler在其有关月球运行理论的著作中，创立了广泛用于求初值问题(1)的数值解的方法，次年又把它推广到二阶方程。欧拉的想法如下：选择步长

h?0，然后在x0?x?x0?h情况下用解函数的切线

y??x?0?x l?x??0代替解函数。这样对于点x1?f,0x?0 y?x0?h就可以得到

. y1?y0?hf?x,0y?0在点?x1,y1?重复如上的程序再次计算新的方向，就会得到所谓的递推公式：

xm?1?xm?h，ym?1?ym?hf?xm,ym?,

这就是Euler方法。由此，再通过连接所有这些切线得到的函数被称为Euler折线。

- 4 -

如果我们令h?0，这些折线就会越来越接近解函数。 (一)、欧拉法

Euler方法是最简单的一步法，它是一阶的，精度较差，但公式很简单，即

yn?1?yn?hy?xn,yn? ?n?0,1,?2? (3)

Euler方法的几何意义在数值计算思想中已经体现出来了，实际上就是用过已知点的折线来近似代替过此点的积分曲线。因此，这种方法又称为折线法。在Euler法中，数值解的误差首先是由差商代替导数引起的，这种近似替代所产生的误差称为截断误差。另外，计算过程中还会由于数值的舍入产生另一种误差——舍入误差。显然只有当初产生的误差在以后各步的计算中不会无限制扩大时，即当初始误差充分小时，以后各步的误差也可以充分小，Euler法才具有实用价值。收敛性、截断误差估计与稳定性闷题是常微分方程各种数值解法研究中必须考虑的基

本问题。显然这些问题在Euler法中是得到验证的，详见下面例子分析。 在xy平面上，微分方程(1)的解y=y(x)称作它的积分曲线。积分曲线上一点?x,y? 线斜率等于函数f(x,y)的值，如果按函数f(x,y)在xy平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致，基于上述几何解释，从初始点p0(x0,y0)出发，先依方向场在该点的方向推进到x?x1上一点p1，然后再从p1依方向场的方向推进到x?x2上一点p2，循此前进推出一条折线

p0p1p2…，一般地，设已做出该折线的顶点pn，过pn(xn,yn)依方向场的方向再推

进到pn?1(xn?1,yn?1)，显然两个顶点pn，pn?1的坐标有关系

yn?1?yn?f?xn,yn?

xn?1?xn即

yn?1?yn?hf?xn,yn? (4) 这就是著名的欧拉(Euler)公式。若初值y0已知，则依公式(4)可逐步算出

y1?y0?hf?x0,y0?,y2?y1?hf?x1,y1?, ?例1 求解初值问题

- 5 -

??y??y?2x/y,(0?x?1)?y(0)?1, (5)

解 欧拉公式的具体形式为

取步长h?0.1，计算结果如下表：

表1 计算结果对比

0.1 0.2 0.3 0.4 0.5 1.1000 1.1918 1.2774 1.3582 1.4351 1.0954 1.1832 1.2649 1.3416 1.4142 0.6 0.7 0.8 0.9 1.0 1.5090 1.5803 1.6498 1.7178 1.7848 1.4832 1.5492 1.6125 1.6733 1.7321 初值问题(2.1.2)有解y?1?2x，按这个解析式子算出的准确值y(xn)同近似值

yn一起列在表1， 两者相比较可以看出欧拉方法的精度很差。

(二)、向后EuIer方法

向后Euler方法和Euler方法差不多，只是把y??xn?1?用 去代替，这时计算公式为

y?xn?1??y?xn?h

??yn?1?y?hf?xn?1,yn?1?, ? ???(n?0,1,2?) （6）

??y?x0??y0,向后Euler方法的总体截断误差也是一阶的，因此向后Euler方法是收敛的。这里需要指出它与Euler方法的一个很大不同之处，Euler方法是显式方法，即yn?1- 6 -