常微分方程数值解法

第七章 常微分方程数值解法常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。

另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。

因此，有必要探讨常微分方程初值问题的数值解法。

本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。

第一节 欧拉法求解常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy（1）的数值解，就是寻求准确解)(x y 在一系列离散节点 <<<<<n x x x x 210 上的近似值 ,,,,,210n y y y y{}n y 称为问题的数值解，数值解所满足的离散方程统称为差分格式，1--=i i ix x h 称为步长，实用中常取定步长。

显然，只有当初值问题(1)的解存在且唯一时，使用数值解法才有意义，这一前提条件由下 面定理保证。

定理 设函数()y x f ,在区域+∞≤≤-∞≤≤y b x a D ,:上连续，且在区域D 内满足李普希兹(Lipschitz)条件，即存在正数L ，使得对于R 内任意两点()1,y x 与()2,y x ,恒有()()2121,,y y L y x f y x f -≤-则初值问题(1)的解()x y 存在并且唯一。

一、欧拉法(欧拉折线法)若将函数)xy 在点nx处的导数()n x y '用两点式代替， 即()hx y x y x y n n n )()(1-≈'+，再用n y 近似地代替()n x y ,则初值问题(1)变为⎩⎨⎧==++=+ ,2,1,0),()(001n x y y y x hf y y n n n n(2)(2)式就是著名的欧拉(Euler)公式。

以上方法称 为欧 拉法或欧拉折线法。

欧拉公式有明显的几何意义。

从几何上看，求解初值问题(1)就是xy 平面上求一条通过点()00,y x 的曲线()x y y =，并使曲线上任意一点()y x ,处的切线斜率为()y x f ,。

欧拉公式的几何意义就是从点)000,y x P 出发作一斜率为()00,y x f 的直线交直线1x x =于点()111,y x P ，1P 点的纵坐标1y 就是()1x y 的近似值；再从点1P 作一斜率为()11,y x f 的直线2x x =交直线于点()222,y x P , 2P 点的纵坐标2y 就是的近似值()2x y ；如此继续进行，得一条折线 210P P P 。

该折线就是解()x y y =的近似图形，如图7-1。

图7-1欧拉法的其它几种解释：1．假设()x y 在n x 附近展开成泰勒级数()()()()()()()()+''++=+''+'+=+n n n n n n n n x y h x y x hf x y x y hx y h x y x y 2,2221取h 的线性部分，并用n y 作为()n x y 的近似值，得 ()n n n n y x hf y y ,1+=+2． 对方程()y x f dxdy,=两边从n x 到1+n x 积分，得()dxx y x f x y x y n nx x n n ⎰+=-+1))(,()(1 (3)用 矩形公式计算上式右侧积分，即 ()()x x x x x d x y x f dx x y x f n nn n⎰⎰++≈11,))(,(并用n y 作为的近似值()n x y ，得()n n n n y x hf y y ,1+=+ 故欧拉法也称为矩形法。

二、改进的欧拉法(梯形法)欧拉法形式简单，低，特别当曲线y=y(x)计算方便，但精度比较的曲率较大时，欧拉法的效果更差。

为了达到较高精度的计算公式，对欧拉法进行改进，用梯形公式计算(3)式 右侧积分，即()()()()[]11,,2))(,(1+++≈⎰+n n n n x x x y x f x y x f h dx x y x f n n并用y n 作为y(x n )的近似值，得到改进 的欧拉公式()()[]111,,2+++++=n n n n n n y x f y x f h y y(4)上述方法称为改进的欧拉法，也称为梯形法。

不难发现，欧拉公式 是 关于y 1+n 的显式，即只要已知y n ,经过一次计算便可得y 1+n 的值，而改进的 欧拉公式是以y 1+n 的隐式方程给出，不能直接得到y 1+n 。

隐式方程(4)通常用 迭代法求解，而迭代过程的实质是逐步显式化。

先用欧拉 公式()()n n n n y x hf y y ,01+=+给出y 1+n 的迭代初值，然后 再用改进的欧拉公式(4)进行迭代，即有[]⎪⎪⎩⎪⎪⎨⎧=++=+=+++++ ,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 nk n n n n n (5) 迭代过程进行到连续两次迭代结果之差的绝对值小于给定的精度ε即()εkn k n y y 111+++-为止，这时取()111+++=k n n y y然后再转入下一步计算。

下面讨论(){}kn y 1+是否收敛；若收敛，它的极限是否满足(4)式。

假设f(x,y)满足李普希兹条件 ()()()2121,,y y L y x f y x f -≤- 则()()()()()()()()()()()011122111211111111111222,,2++-+-+-++-+++++++-⎪⎭⎫⎝⎛≤≤-⎪⎭⎫ ⎝⎛≤-≤-=-n n k n k n k n k n k n n k n n k n k n y y hL y y hL y y hL y x f y x f h y y由此可见，只要12<hL (这里只要步长h 足够小即可)，当k →∞时，有02→⎪⎭⎫⎝⎛khL ，所以(){}kn y 1+收敛。

又因为f(x ,y)对y 连续，当k →∞时，对等式()()()()[]k n n n nn k n y x f y xf h y y 1111,,2++++++=两端取极限， 得()()[]111,,2+++++=n n n n n n y x f y x f h y y因 此，只要步长h 足够小，就可保证(){}kn y 1+收敛到满足(4)式的1+n y。

三、预估一校正法改进的欧拉公式在实际计算时要进行多次 迭代，因而计算量较大。

在实用上，对于改进的欧拉公式(5)只迭代一次，即先用欧拉公式算出y 1+n 的预估值y)0(1+n ，再用改进的欧拉公式(4)进行一次迭代得到 校正值y 1+n ，即()()[]⎪⎩⎪⎨⎧=++=+=++++ ,2,1,0,,,2),(111)0(1n y x f y x f h y y y x hf y y n n n n n n n n n n (6)预估—校正公式也常写成下列形式：(),2,1,0,,),(2121121211=⎪⎪⎩⎪⎪⎨⎧++==++=+n k y h x hf k y x hf k k k y y n n n n n n (7)四、公式的截断误差定义 若某种微分方程数值解 公式的截断误差是O(h1+k )，则称这种方法是k 阶方法。

为了简化分析，在进行误差分析时，我们假设前一步的结果是准确的，即在 y n =y(x n )的前提下，考虑用y 1+n 作为y(x 1+n )的近似值而产生的截断误差，这种误差称为局部截断误差。

由泰勒公式()()()()()+''+'+=+=+n n n n n x y hx y b x y h x y x y !221对于欧拉公式，有()()()n n n n n n x y h x y y x hf y y '+=+=+,1 于是()()211hO y x y n n =-++则欧拉公式的截断误差为O(2h ),所以 欧拉法是一阶方法。

对于预估—校正公式，有()()()()()()()()()()()[]()()()()()()()[]+'++=+++=++=++='==n n y n n nxn n n n y n n xn n n n n n n n n x y x f x y x y xf hx y x hf x y x f k x y x hf x y x f h k x y h x hf k y h x hf k x y h y x hf k ,,,,,,,,,211121而()()()()()()()()()x y x f x y x y x f x y x y x f x y y x ,,,⋅'+=''='于是()()+''+'=n n x y h x y h k 22因此()()()+''+'+=++=+n n n n n x y hx y h x y k k y y 221212211所以y(x 1+n )- y 1+n = O(3h )则预估—校正公式的截断误差为O(3h ),也即预估—校正法是二阶方法。

可以证明，改进的欧拉公式与预估—校正公式的截断误差相同，均为O(3h )。

这里略去证明。

例 1：在区间［0,1］上以h=0.1为步长，分别用欧拉法与预估— 校正法求初值问题⎪⎩⎪⎨⎧=-=1)0(2y y x y dx dy的数值解解：该方程为贝努利方程，其精确解为x y 21+=欧拉公式的具体形式为⎪⎪⎭⎫ ⎝⎛-+=+nnnn n y x y h y y 21预估—校正公式的具体形式为()⎪⎪⎪⎩⎪⎪⎪⎨⎧⎪⎪⎭⎫⎝⎛++-+=-=++=+11212112)2(2121k y h x k y hf k y x y hf k k k y y n n nn nn n n 取步长1,0,1.000===y x h ,计算结果见下表：表7-1表7-1近似解与准确解比较，欧拉法的结果大致只有两位有效数字，而预估—校正法的结果则 有3位有效数字。

第二节 龙格—库塔法前面讨论的欧拉法与改进的欧拉法都是一步 法，即计算y 1+n 时，只用到前一步值。

龙格—库塔(Runge-Kut ta)法(简称为 R-K 方法)是一类高精度的一步法，这类方法与泰勒级数法有着密 切的关系。

一、泰勒级数法设有初值问题⎪⎩⎪⎨⎧==00)(),(y x y y x f dx dy由 泰勒展开式()()()()()()()()121!!2+++++''+'+=+=k n k kn n n n n hO x yk hx y hx y h x y h x y x y （1）若令()()()()n k kn n n x yk hhx y h x y y !!221+++'+=+ (2)则()()111+++=-k n n hO y x y即公式(2)为k 阶方法。