数值分析常微分方程数值解

许多实际问题的数学模型是微分方程或微分方程的定解问题。

如物体运动、电路振荡、化学反映及生物群体的变化等。

常微分方程可分为线性、非线性、高阶方程与方程组等类；线性方程包含于非线性类中，高阶方程可化为一阶方程组。

若方程组中的所有未知量视作一个向量，则方程组可写成向量形式的单个方程。

因此研究一阶微分方程的初值问题⎪⎩⎪⎨⎧=≤≤=0)(),(y a y bx a y x f dxdy， (9-1) 的数值解法具有典型性。

常微分方程的解能用初等函数、特殊函数或它们的级数与积分表达的很少。

用解析方法只能求出线性常系数等特殊类型的方程的解。

对非线性方程来说，解析方法一般是无能为力的，即使某些解具有解析表达式，这个表达式也可能非常复杂而不便计算。

因此研究微分方程的数值解法是非常必要的。

只有保证问题(9-1)的解存在唯一的前提下，研究其数值解法或者说寻求其数值解才有意义。

由常微分方程的理论知，如果(9-1)中的),(y x f 满足条件（1）),(y x f 在区域} ),({+∞<<∞-≤≤=y b x a y x D ，上连续； （2）),(y x f 在上关于满足Lipschitz 条件，即存在常数，使得y y L y x f y x f -≤-),(),(则初值问题(9-1)在区间],[b a 上存在惟一的连续解)(x y y =。

在下面的讨论中，我们总假定方程满足以上两个条件。

所谓数值解法，就是求问题(9-1)的解)(x y y =在若干点b x x x x a N =<<<<= 210处的近似值),,2,1(N n y n =的方法。

),,2,1(N n y n =称为问题(9-1)的数值解，n n x x h -=+1称为由到1+n x 的步长。

今后如无特别说明，我们总假定步长为常量。

建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： (1) 用差商近似导数在问题(9-1)中，若用向前差商hx y x y n n )()(1-+代替)(n x y '，则得)1,,1,0( ))(,()()(1-=≈-+N n x y x f hx y x y n n n n n)(n x y 用其近似值代替，所得结果作为)(1+n x y 的近似值，记为1+n y ，则有 1(,) (0,1,,1)n n n n y y hf x y n N +=+=-这样，问题(9-1)的近似解可通过求解下述问题100(,) (0,1,,1)()n n n n y y hf x y n N y y x +=+=-⎧⎨=⎩(9-2)得到，按式(9-2)由初值经过步迭代，可逐次算出N y y y ,,21。

此方程称为差分方程。

需要说明的是，用不同的差商近似导数，将得到不同的计算公式。

(2) 用数值积分法将问题(9-1)中的微分方程在区间],[1+n n x x 上两边积分，可得)1,,1,0( ))(,()()(11-==-⎰++N n dx x y x f x y x y n nx x n n (9-3)用1+n y ,分别代替)(1+n x y ,)(n x y ，若对右端积分采用取左端点的矩形公式，即),())(,(1n n x x y x hf dx x y x f n n≈⎰+同样可得出显式公式(9-2)。

类似地，对右端积分采用其它数值积分方法，又可得到不同的计算公式。

(3) 用Taylor 多项式近似。

把1()n y x +在点处Taylor 展开，取一次多项式近似，则得2121()()()()2!()(,())() [,]2!n n n n n n n n h y x y x hy x y hy x hf x y x y x x ξξξ++'''=++''=++∈ 设1h ，略去余项，并以代替()n y x ，便得 1(,)n n n n y y hf x y +=+以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式。

其中Taylor 展开法，不仅可以得到求数值解的公式，而且容易估计截断误差。

上面我们给出了求解初值问题(9-1)的一种最简单的数值公式(9-2)。

虽然它的精度比较低，实践中很少采用，但它的导出过程能较清楚地说明构造数值解公式的基本思想，且几何意义明确，因此它在理论上仍占有一定的地位。

1简单的数值方法和基本概念1.1 Euler 法与向后Euler 法一、Euler 法Euler 方法就是用差分方程初值问题10(,) (0,1,,1)()n n n n y y hf x y n N y y a +=+=-⎧⎨=⎩(9-4)的解来近似微分方程初值问题(9-1)的解，即由公式(9-4)依次算出()n y x 的近似值(1,2,)n y n =。

从几何上看，微分方程(,)y f x y '=在xoy 平面上确定了一个向量场：点(,)x y 处的方向斜率为(,)f x y 。

问题(9-1)的解()y y x =代表一条过点00(,)x y 的曲线，称为积分曲线，且此曲线上每点的切向都与向量场在这点的方向一致。

从点000(,)P x y 出发，以00(,)f x y 为斜率作一直线段，与直线1x x =交于点111(,)P x y ，显然有1000(,)y y hf x y =+，再从出发，以11(,)f x y 为斜率作直线段推进到2x x =上一点222(,)P x y ，其余类推，这样得到解曲线的一条近似曲线，它就是折线012P PP 。

因此Euler 方法又称为Euler 折线法。

二、向后Euler 法在微分方程离散化时，用向后差商代替导数，即11()()()n n n y x y x y x h++-'≈，则得到如下差分方程11100(,) (0,1,,1)()n n n n y y hf x y n N y y x +++=+=-⎧⎨=⎩(9-5)用这组公式求问题(9-1)的数值解称为向后Euler 法。

向后Euler 法与Euler 法形式上相似，但实际计算时却复杂得多。

Euler 法计算1+n y 的公式中不含有1+n y ，这样的公式称为显式公式；向后Euler 法计算1+n y 的公式中含有1+n y ，称为隐式公式。

显式公式与隐式公式各有特点。

显式公式的优点是使用方便，计算简单，效率高。

其缺点是计算精度低，稳定性差；隐式公式正好与它相反，它具有计算精度高，稳定性好等优点，但求解过程很复杂，一般采用迭代法。

为了结合各自的优点，通常将显式公式与隐式公式配合使用，由显式公式提供迭代初值，再经隐式公式迭代校正。

上面隐式公式中，在求解1n y +时，1,n n y x +为已知，1n y +是方程111(,)n n n n y y hf x y +++=+的根。

一般说来，这是一个非线性方程，因此我们通过构造简单迭代法来求解。

迭代格式为(0)1(1)()111(,)(,) (0,1,2,)n n n n k k n n n n y y hf x y y y hf x y k +++++⎧=+⎪⎨=+=⎪⎩ 由于(,)f x y 满足Lipschitz 条件，所以(1)()()11111111(,)(,)k k k n n n n n n n n yy h f x y f x y hL yy +++++++++-=-≤-由此可知，只要1hL <，迭代法就收敛到解1n y +。

1.2 梯形公式利用数值积分方法将微分方程离散化时，若用梯形公式计算式(9-3)中右端积分，即111(,())[(,())(,())]2n nx n n n n x hf x y x dx f x y x f x y x +++≈+⎰ 并用1,n n y y +代替1(),()n n y x y x +，则得计算公式111[(,)(,)]2n n n n n n hy y f x y f x y +++=++(9-6)这就是求解初值问题(9-1)的梯形公式。

梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为(0)1(1)()111(,)(,)(,) (0,1,)2n n n n k k n n n n n n y y hf x y h y y f x y f x y k +++++⎧=+⎪⎨⎡⎤=++=⎪⎣⎦⎩ 由于函数(,)f x y 关于满足Lipschitz 条件，所以(1)()()11111111(,)(,)22k k k n n n n n n n n h hL y y f x y f x y y y +++++++++-=-≤- 其中为Lipschitz 常数。

因此，当012hL<<时，迭代法就收敛到解1n y +。

1.3 局部截断误差与方法的精度为了刻画近似解的准确程度，引入局部截断误差与方法精度的概念。

定义 假设在某一步的近似解是准确的，即()n n y y x =（这个假设称为局部化假设）。

在此前提下，用某公式推算所得1n y +，我们称111()n n n R y x y +++=-为该公式（即该方法）的局部截断误差。

定义9.2 如果某种方法的局部截断误差是1111()()p n n n R y x y O h++++=-=则称该方法是阶方法，或具有阶精度。

显然越大，方法的精度越高。

1）Euler 法的截断误差假设问题的解)(x y 充分光滑，且前步计算结果是精确的，即)(i i x y y =，()(,())i i i y x f x y x '=)(n i ≤于是Euler 法的截断误差是11111()()(,)()()()n n n n n n n n n n R y x y y x y hf x y y x y x hy x +++++=-=--'=--23()()2n h y x O h ''=+(9-7) 这里2()2n h y x ''称为局部截断误差主项。

显然21()n R O h +=。

2）向后Euler 法的截断误差。

计算公式是),(111++++=n n n n y x hf y y将(,)f x y 对用微分中值定理，有1111111(,)(,())(,)(())n n n n y n n n f x y f x y x f x y y x η+++++++=+-（在1n y +与1()n y x +之间）将11(,())n n f x y x ++在处Taylor 展开2111(,())()()()()n n n n n f x y x y x y x hy x O h +++''''==++于是231111()()()(,)(())()n n n n y n n n y y x hy x h y x hf x y y x O h η++++'''=+++-+将方程的解作Taylor 展开231()()()()()2n n n n h y x y x hy x y x O h +'''=+++因此2311111()()(,)(())()22n n n y n n n h hy x y y x f x y y x O h η+++++''-=---+故2311112311()()()1(,)21(,)()()2n n n n y n y n n h R y x y y x O h hf x h hf x y x O h ηη+++++⎡⎤''=-=-+⎢⎥-⎣⎦⎡⎤''⎡⎤=++-+⎢⎥⎣⎦⎣⎦23()()2n h y x O h ''=-+(9-8)3）梯形法的计算公式是111[(,)(,)]2n n n n n n hy y f x y f x y +++=++将(,)f x y 对用微分中值定理，有1111111(,)(,())(,)(())n n n n y n n n f x y f x y x f x y y x η+++++++=+-（在1n y +与1()n y x +之间）将11(,())n n f x y x ++在处Taylor 展开23111(,())()()()()()2n n n n n n h f x y x y x y x hy x y x O h +++'''''''==+++图9-3于是2341111()()()()(,)(())()242n n n n n y n n n h h hy y x hy x y x y x f x y y x O h η++++''''''=++++-+将方程的解作Taylor 展开2341()()()()()()26n n n n n h h y x y x hy x y x y x O h +''''''=++++因此3411111()()(,)(())()122n n n y n n n h hy x y y x f x y y x O h η+++++'''-=---+ 故341111341341()()()1(,)12[1(,)]()()12()()12n n n ny n y n n n h R y x y y x O h hf x h hf x y x O h h y x O h ηη+++++⎡⎤'''=-=-+⎢⎥-⎣⎦⎡⎤'''=++-+⎢⎥⎣⎦'''=-+(9-9)所以梯形法是二阶方法。