5 常微分方程组与高阶常微分方程的数值解法

§5 常微分方程组与 高阶常微分方程的数值解法

为了简单起见，在此仅研究含有两个一阶常微分方程的初值问题

111210102

2122020(,,),(),(,,),()y f x y y y x y y f x y y y x y '==⎧⎨

'==⎩ （5.1） 若记

12y y ⎛⎫

= ⎪⎝⎭y

12(,)(,)(,)f x x f x ⎛⎫

= ⎪

⎝⎭y f y y

10020y y ⎛⎫= ⎪

⎝⎭y

10020()()()y x x y x ⎛⎫= ⎪⎝⎭

y

则（5.1）可以改写为向量的形式

00(,)()x x '=⎧⎨

=⎩y f y y y （5.2）

相应地，将前述求解常微分方程初值问题的公式写成向量的形式，就得到常微分方程组的数值解法。 例如，求解（5.1）的四阶经典R-K 方法为

1,11131112142,1223212224226k k k k y y K K K K h y y K K K K ++⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫

⎛⎫=++++ ⎪⎢⎥ ⎪ ⎪ ⎪ ⎪ ⎪

⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦

其中

1121121221(,,)(,,)k k k k k k f x y y K f x y y K ⎛⎫

⎛⎫=

⎪ ⎪⎝⎭⎝⎭

111122112222111221(,,)222

(,,)222k k k k k k h h h f x y K y K K K h h h f x y K y K ⎛⎫+++ ⎪⎛⎫=

⎪ ⎪ ⎪⎝⎭+++ ⎪⎝⎭

111222213232112222(,,)222

(,,)222k k k k k k h h h f x y K y K K K h h h f x y K y K ⎛⎫+++ ⎪⎛⎫=

⎪ ⎪ ⎪⎝⎭+++ ⎪⎝⎭

111322314211322324(,,)(,,)k k k k k k f x h y hK y hK K f x h y hK y hK K +++⎛⎫

⎛⎫=

⎪ ⎪+++⎝⎭⎝⎭

注：教材中1122(),()k k k k y y x y y x ==不正确，它们仅当0k =时成立（即初始条件）。

对于高阶微分方程，可以化成微分方程组后再求解。例如二阶常微分方程初值问题

0000(,,),(),()y f x y y y x y y x y '''=⎧⎪

=⎨⎪''=⎩

（5.3）

令

12,y y y y '==，

则（5.3）可以化为

12100102

1220020,()(,,),()()y y y x y y y f x y y y x y x y '===⎧⎨

''===⎩

§6 边值问题的数值解法

本节用最简单的二阶常微分方程问题为例，说明边值问题的数值解法。

(,,),[,],(),().y f x y y x a b y a y b αβ'''=∈⎧⎪

=⎨⎪=⎩

用中心差商代替导数，上述问题可以转化为方程组，然后求解。这种求解方法称为差分法。过程如下： 记

00,(1,2,,),k b a h x x kh k n x a

n

-==+== ， 将一阶和二阶中心差商公式

2

11()()()()

2k k k y x y x y x O h h

+--'=+ 2

112()2()()()()k k k k y x y x y x y x O h h

-+-+''=+ 中的余项2

()O h 略去，并用11,,k k k y y y -+分别替代

11(),(),()k k k y x y x y x -+，代入方程，便得

2

111102(,,),21,2,,1,,

.

k k k k k k k n y y y y y h f x y h k n y y αβ+-+--⎧-+=⎪⎪⎪=-⎨

⎪=⎪

=⎪⎩ 注意到0,n y y 为已知，这是一个含有1n -个未知数，

1n -个方程的方程组。

一般情况下，上述方程组是非线性的，求解比较困难。如果(,,)f x y y '是y 与y '的线性函数，则方程

组为线性方程组。

例2 求边值问题

,[0,1],(0)0,

(1)1

y y x x y y ''-=∈⎧⎪

=⎨⎪=⎩

取步长0.1h =，求数值解。

解 本例中，

0(,,)0,1,0,1110

0.1,0,1,,n k f x y y y x a b y y n h x k k n

'=+======== 分别取1,2,,1k

n =- ，得

2

21012

3212210

98920.1(0.1)20.1(0.2)20.1(0.9)

y y y y y y y y y y y y ⎧-+=+⎪-+=+⎪⎨⎪⎪-+=+⎩

将0100,1y y ==代入上式，整理后得

21321982.010.0012.010.0022.010.0091

y y y y y y y -=⎧⎪-+=⎪⎨⎪⎪-+=-⎩ 上述方程组为对角占优的三对角方程组，用追赶法求

解，得