081数值计算方法—常微分方程组

科学计算—理论、方法

及其基于MATLAB 的程序实现与分析

微分方程（组）数值解法

§1 常微分方程初值问题的数值解法

微分方程（组）是科学研究和工程应用中最常用的数学模型之一。如揭示质点运动规律的Newton 第二定律：

()()()⎪⎩⎪⎨⎧'='==0

00022x t x x t x t F dt x

d m （１） 和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题：

()

()0

0y y t f ay dt

dy

=+= （２） 的解为：

()()()τττd f e y e t y t

t a at ⎰-+=00 （３）

但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：

１、绝大多数的常微分方程和常微分方程组得不到（有限形式的）解析解；

２、实际应用中往往只需要知道常微分方程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）；

如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。

一般的一阶常微分方程（组）的初值问题是指如下的一阶常微分方程（组）的定解问题：

()()0

00,y t y t t t y t F dt

dy

f

=≤≤= （７）

其中

()()()()⎪⎪⎪⎪

⎪⎭

⎫ ⎝⎛=t y t y t y t y n 21 （８）

()()()()⎪⎪⎪⎪

⎪⎭

⎫ ⎝⎛=y t f y t f y t f y t F n ,,,,21

（９） 常微分方程（组）的初值问题通常是对一动态过程（动态系统、动力系统）演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。

§１.1 常微分方程（组）的Cauch 问题数值解法概论

假设要求在点（时刻）k t ，n k ,,2,1 =处初值问题（７）的解

()00,,y t t y y =的（近似）值,如果已求得k t 时刻的值()k t y 或它的近似值k y (如0t 时刻的值()0t y )，那么将式（７）的两端在区间[]1,+k k t t 上积分

()()()()dt t y t f t y t y dt dt

dy

k k

k k

t t k k t t ⎰⎰

++=-=+11

,1 (10)

可得

()()()()dt t y t f t y t y k k

t t k k ⎰

++

=+1

,1 (11)

或

()()()dt t y t f y y t y k k

t t k k k ⎰

++

=≈++1

,11 (12)

显然,为了利用式（11）或（12）求得()1+k t y 的精确值（近似值1+k y ），必须计算右端的积分，这是问题的关键也是难点所在，如前所述，一般得不到精确的公式解，因此需要采用数值积分的方法求其近似解, 可以说,不同的式值积分方法将给出不同的Cauch 问题的数值解法。 §1.2 最简单的数值解法——Euler 方法

假设要求在点（时刻）kh t t k +=0，n

t t h f 0

-=

，n k ,,2,1 =处初值问

题（７）的解()t y y =的近似值。首先对式（７）的两端积分，得

()()()()dt t y t f t y t y dt dt

dy

k k k k

t t k k t t ⎰⎰

++=-=+11

,1 （13） 对于式（13）的右边，如果用积分下限k t t =处的函数值()()k k t y t f ,代替被积函数作积分（从几何上的角度看，是用矩形面积代替曲边梯形面 积），则有

()()()()()[]k k t t

k k t y t hf dt t y t f t y t y k k

,,1

1≈=-⎰++ （14）

进而得到下式给出的递推算法—Euler 方法

()()0

01,,2,1,y t y n

k y t hf y y k k k k ==+=+ （15）

例１ 用Euler 方法解如下初值问题，取3.0=h ，

()1

03

0sin 23=≤≤-=y t t y y dt

dy

解：由（15）得

()1

010

,,2,1sin 6.03.13

1==-=+y k t y y y k

k k k

结果如下：

如果取1.0=h ，其结果如下图所示：

Euler_Method

§1.3 改进的Euler 方法

对于（15）的右边，如果被积函数用积分限k t t =和1+=k t t 处的函数值的算术平均值代替（几何上，是用梯形面积代替曲边梯形面积），则有

()()()()()()()()[]111,,2

,1++++≈

=-⎰

+k k k k t t k k t y t f t y t f h

dt t y t f t y t y k k

（16） 进而得到下式给出的递推算法：

()()[]()0

0111,,2,1,,2

y t y n

k y t f y t f h

y y k k k k k k ==++=+++ （17）