第七章 常微分方程数值解法
常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。因此，有必要探讨常微分方程初值问题的数值解法。本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。
第一节 欧拉法
求解常微分方程初值问题
?dy??f(x,y) ?dx??y(x0)?y0 （1）
的数值解，就是寻求准确解y(x)在一系列离散节点x0?x1?x2???xn?? 上的近似值 y0,y1,y2,?,yn,?
?yn?称为问题的数值解，数值解所满足的离散方程统称为差分格式，hi称为步长，实用中常取定步长。
?xi?xi?1显然，只有当初值问题(1)的解存在且唯一时，使用数值解法才有意义，这一前提条件由下 面定理保证。
定理 设函数f?x,y?在区域D:a?x?b,???y???
上连续，且在区域D内满足李普希兹(Lipschitz)条件，即存在正数L，使得对于R内任意两点?x,y1?与?x,y2?,恒有的解y?x?存在并且唯一。 一、欧拉法(欧拉折线法)
若将函数
yxf?x,y1??f?x,y2??Ly1?y2则初值问题(1)
?在点xhn?处的导数y?xn?用两点式代替， 即
y??xn??y(xn?1)?y(xn)，再用yn近似地代替y?xn?,则初值问题(1)变为
?yn?1?yn?hf(xn?yn) (2) ?y?y(x),n?0,1,2,?0?0
1
(2)式就是著名的欧拉(Euler)公式。以上方法称 为欧 拉法或欧拉折线法。?欧拉公式有明显的几何意义。从几何上看，求解初值问题(1)就是xy平面上求一条通过点?x0,y0?的曲线y?y?x?，并使曲线上任意一点?x,y?处的切线斜率为
f?x,y?。 欧拉公式的几何意义就是从点
P0x0,y0?出发作一斜率为f?x0,y0?的
直线交直线x?x1于点P1?x1,y1?，P1点的纵坐标y1就是y?x1?的近似值；再从点P1作一斜率为f?x1,y1?的直线x?x2交直线于点P2?x2,y2?, P2点的纵坐标y2就是的近似值y?x2?；如此继续进行，得一条折线P0P1P2?。该折线就是解y?y?x?的近似图形，如图7-1。
图7-1
欧拉法的其它几种解释：
1．
假设y?x?在xn附近展开成泰勒级数
y?xn?1??y?xn??hy??xn???y?xn??hf?xn,y?xn???hh222y???xn???2y???xn???
取h的线性部分，并用yn作为y?xn?的近似值，得 yn?1?yn?hf?xn,yn?
dy2． 对方程
dx?f?x,y?两边从xn到xn?1积分，得
y(xn?1)?y?xn???xn?1xnf(x,y(x))dx (3)
2
用 矩形公式计算上式右侧积分，即 ?xn?1xnf(x,y(x))dx??xn?1xnf?x,y?x??dx
并用yn作为的近似值y?xn?，得yn?1?yn?hf?xn,yn? 故欧拉法也称为矩形法。 二、改进的欧拉法(梯形法)
欧拉法形式简单，低，特别当曲线y=y(x)计算方便，但精度比较的曲率较大时，欧拉法的效果更差。为了达到较高精度的计算公式，对欧拉法进行改进，用梯形公式计算(3)式 右侧积分，即
?xn?1xnf(x,y(x))dx?h2?f?xn,y?xn???f?xn?1,y?xn?1???
并用yn作为y(xn)的近似值，得到改进 的欧拉公式
yn?1?yn?h2?f?xn,yn??f?xn?1,yn?1?? (4)
上述方法称为改进的欧拉法，也称为梯形法。
不难发现，欧拉公式 是 关于yn?1的显式，即只要已知yn,经过一次计算便可得yn?1的值，而改进的 欧拉公式是以yn?1的隐式方程给出，不能直接得到yn?1。隐式方程(4)通常用 迭代法求解，而迭代过程的实质是逐步显式化。
先用欧拉 公式
yn?1?yn?hf?xn,yn??0?
给出yn?1的迭代初值，然后 再用改进的欧拉公式(4)进行迭代，即有
(0)?yn?1?yn?hf(xn,yn)?h?(k?1)(k)y?y?f(xn,yn)?f(xn?1,yn?1)?n?1n (5) 2?k?0,1,2,?????迭代过程进行到连续两次迭代结果之差的绝对值小于给定的精度?即
yn?1?yn?1???k?1?k
为止，这时取
yn?1?yn?1?k?1?
然后再转入下一步计算。
3
???y?下面讨论是否收敛；若收敛，它的极限是否满足(4)式。
kn?1假设f(x,y)满足李普希兹条件 f?x,y1??f?x,y2??L?y1?y2? 则
yn?1?yn?1???k?1??k?h2f?xn?1,yn?1??f?xn?1,yn?1?k??k?1??hL2yn?1?yn?12?k??k?1??hL?????2?yn?1?yn?12?k?1?k?2?hL??????2??yn?1?yn?1?1??0?
?hL??k→∞时，有??2?k由此可见，只要
kn?1hL2?1(这里只要步长h足够小即可)，当
?0，
???y?所以收敛。又因为f(x ,y)对y连续，当k→∞时，对等式
yn?1?k?1??yn?h2?f?xn,yn??fxn?1,yn?1??k???
两端取极限， 得
yn?1?yn?h2?f?xn,yn??f?xn?1,yn?1??kn?1
n?1???y?因 此，只要步长h足够小，就可保证收敛到满足(4)式的y。
三、预估一校正法
改进的欧拉公式在实际计算时要进行多次 迭代，因而计算量较大。在实用上，对于改进的欧拉公式(5)只迭代一次，即先用欧拉公式算出yn?1的预估值y
(0)n?1，再用改进的欧拉公式(4)进行一次迭代得到 校正值yn?1，即
(0)?yn?1?yn?hf(xn,yn)?h?yn?1?yn??f?xn,yn??f?xn?1,yn?1??,n?0,1,2,??2? (6)
预估—校正公式也常写成下列形式：
4
11?y?y?k?k2n1?n?122?k1?hf(xn,yn),n?0,1,2,???k?hf?x?h,y?k?nn1?2? (7)
四、公式的截断误差
定义 若某种微分方程数值解 公式的截断误差是O(h
k?1)，则称这种方法是
k阶方法。为了简化分析，在进行误差分析时，我们假设前一步的结果是准确的，即在 yn=y(xn)的前提下，考虑用yn?1作为y(xn?1)的近似值而产生的截断误差，这种误差称为局部截断误差。由泰勒公式
y?xn?1??y?xn?h??y?xn??by??xn??h22!y???xn???
对于欧拉公式，有
yn?1?yn?hf?xn,yn??y?xn??hy??xn? 于是
y?xn?1??yn?1?Oh2??
2则欧拉公式的截断误差为O(h),所以 欧拉法是一阶方法。
对于预估—校正公式，有
k1?hf?xn,yn??hy??xn?k2?hf?xn?h,yn?k1??hf?xn?h,y?xn??k1??hf?xn,y?xn???hf?hf?xn,y?xn???h2?x?xn,y?xn???k1fy?xn,y?xn?????xn?f?x,y?xn???y??xn?fy?xn,y?xn?????
而
y??x??f?x,y?x??y???x??fx?x,y?x???y??x??fy?x,y?x?? 于是
k2?hy??xn??hy???xn???2
因此
yn?1?yn?12k1?312k2?y?xn??hy??xn??h22y???xn???
3所以y(xn?1)- yn?1= O(h)则预估—校正公式的截断误差为O(h),也即预估—
5