(整理)常微分方程数值解法

i.常微分方程初值问题数值解法

i.1 常微分方程差分法

考虑常微分方程初值问题：求函数()u t 满足

(,), 0du f t u t T dt

=<≤ （i.1a ） 0(0)u u = (i.1b)

其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得

121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-∀∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。

通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法--差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点：

0110N N t t t t T -=<<<<=L （i.3） 其中n t hn =，0h >称为步长。

在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商

10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++

如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到

1000(,)u u hf t u -=

一般地，我们有

1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-L 方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t L 上的差分解1,,N u u L 。

下面我们用数值积分法重新导出 Euler 法以及其它几种方法。为此，在区间1[,]n n t t +上积分常微分方程（i.1a ），得

1

1()()(,())n n t n n t u t u t f t u t dt ++=+⎰ （i.5）

用各种数值积分公式计算（i.5）中的积分，便导致各种不同的差分法。例如，若用左矩形公式就得到 Euler 法（i.4）。如果用右矩形公式，便得到下面的：

111Euler (,), 0,1,,1n n n n u u h f t u n N +++=+=-L 隐式方法： （i.6） 类似地，如果用梯形公式，就得到

11111Euler [(,)(,)], 0,1,,12

n n n n n n h u u f t u f t u n N +++++=++=-L 改进的方法 (i.7) 当(,)f t u 关于u 是非线性函数的时候，不能由（i.6）和(i.7)从n u 直接算出1n u +，称这一类方法为隐式，通常采用某种迭代法求解。例如，将一般的隐式方法写成

11(,,)n n n n u F t u u ++= (i.8) 则可以利用如下的迭代法由n u 算出1n u +：

11101(,,), 0,1, k k n n n n n n u F t u u k u u ++++⎧==⎨=⎩

L (i.9) 关于k 的迭代通常只需进行很少几步就可以满足精度要求了。

为了避免对隐式方法进行迭代的麻烦，比如说对于改进的Euler 方法（i .7），可以采用某种预估法近似算出11(,)n n f t u ++，然后再用（i .7）作校正，这就导致所谓预估校正法。下面给出一个例子：

1111111(,)2(,)()2

n n n n n n n n n n n n n u f t u u u hu u f t u h u u u u +-+++++'⎧=⎧⎪⎪'=+⎨⎪⎪⎪'=⎨⎩⎪⎪''=++⎪⎩预估： 校正： （i.10） 这是一个多步法，即计算节点1n t +上的近似值1n u +时，除了用到前一点的近似值n u 之外，还要用到1n u -，甚至可能用到2,n u -L 。而用前面的各种Euler 法计算节点1n t +上的近似值1n u +时，只用到n u ，因此称之为单步法。

下面给出另一个多步法的例子。在区间2[,]n n t t +上积分（i.1a ），得

2

2()()(,())n n t n n t u t u t f t u t dt ++=+⎰

用Simpson 公式(即把被积函数看作二次函数)近似计算积分，便得到

212Milne (4), 0,1,,23

n n n n n h u u f f f n N +++=+++=-L 法： （i.11） 用多步法（i.10）或（i.11）计算时，必须先用某种单步法由0u 计算出1u ，称为造表头。然后再逐次算出23,,,N u u u L 。

一般说来，多步法比Euler 法等简单的单步法精度要高一些。下面我们讨论一类所谓Ronge-Kutta 法。他们是单步法，但是其精度可以与多步法比美。最常用的是下面的标准Ronge-Kutta 法和Gear 法：

112234311234(,)(,)22Lunge-Kutta (,)22(,)(22)6n n n n n n n n n K f t u K h K f t u K h K f t u K f t h u K h u K K K K +=⎧⎪⎪=++⎪⎪⎪=++⎨⎪⎪=++⎪⎪=+++⎪⎩

标准法： （i.12）

(i.13) 从几何上， Ronge-Kutta 法可以粗略地解释为：在区间1[,]n n t t +中选取若干个点(可以重复)1211n k k n t t ττττ-+=<≤≤≤=L ，仅仅利用在区间1[,]n n t t +内可以得到的所有信息，依次给出函数(,())f t u t 在这些点上尽可能精确的的近似值12,,,k K K K L ，然后把它们组合起来，尽可能精确地近似计算（i.5）中的积分。

Ronge-Kutta 法一般的构造方式如下。选定常数k ，令

1(,)n n K f t u =

22211(,)n n K f t h u K αβ=++