常微分方程数值解法

本文格式为Word版，下载可任意编辑——常微分方程数值解法第七章常微分方程数值解法

常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。另外，有些初值问题虽然有初等解，但由于形式太繁杂不便于应用。因此，有必要探讨常微分方程初值问题的数值解法。本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。

第一节欧拉法

求解常微分方程初值问题

?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