常微分方程数值解法

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 -=<<<<= （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 +=+=-方法：（i.4）
从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。

下面我们用数值积分法重新导出 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 +++=+=-隐式方法：（i.6）
类似地，如果用梯形公式，就得到 11111Euler [(,)(,)], 0,1,,12n n n n n n h u u f t u f t u n N +++++=++=-改进的方法 (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 ++++⎧==⎨=⎩ (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 -。

而用前面的各种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,,23n n n n n h u u f f f n N +++=+++=-法： （i.11）
用多步法（i.10）或（i.11）计算时，必须先用某种单步法由0u 计算出1u ，称为造表头。

然后再逐次算出23,,,N u u u 。

一般说来，多步法比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 ττττ-+=<≤≤≤=，仅仅利用在区间1[,]n n t t +内可以得到的所有信息，依
次给出函数(,())f t u t 在这些点上尽可能精确的的近似值12,,
,k K K K ，然后把它们组合起来，尽可能精确地近似计算（i.5）中的积分。

Ronge-Kutta 法一般的构造方式如下。

选定常数k ，令
1(,)n n K f t u =
22211(,)n n K f t h u K αβ=++
33311322(,)n n K f t h u K K αββ=+++
1122,11(,)k n k n k k k k k K f t h u K K K αβββ--=+++++ 11122()n n k k u u h K K K ωωω+=++++ （i.14） 选取这些待定常数,,i ij m αβω的原则是：将（i.14）在(,)n n t u 作Taylor 展开，然后按照h 的幂重新整理，使得
231123112!31n n u u h h h γγγ+=+++
+ （i.15）
与微分方程（i.1a ）的解()u u t =在n t t =处的Taylor 展式
23111()()2!3!
n n n n n u t u t f h f h f h +'''=+++
+ （i.16） 有尽可能多的项重合，即要求 123, , ,
n n n f f f γγγ'''==
= 这里(,)(,)
(,)(,), n n n n n n t u t u df t u f f t u f dt ='==等等。

按照（i.14）构造出的都是显式Runge-Kutta 方法,在每一个i K 的表达式中只出现n u 。

如果允许在某一个i K 的表达式中出现1n u +，则可以导出隐式Runge-Kutta 方法。

i.2 常微分方程组与高阶常微分方程
先来考虑下面的常微分方程组初值问题
1
1122212121100
(,,,,)(,,,,)(,,,,)(0),,(0)m m m
m m m m du f t u u u dt du f t u u u dt du f t u u u dt
u u u u ⎧=⎪⎪⎪=⎪⎪⎨⎪⎪=⎪⎪⎪==⎩ (i.17) 利用向量记号,上式可以改写为
()(,)(0)t t u '=⎧⎨=⎩u f u u (i.18) 上节中各方法都可以直接应用到常微分方程组(i.18)。

例如，Euler 方法成为
1(,), 0,1,
,1n n n n h t n N +=+=-u u f u
再来考虑高阶常微分方程 111211
(,,,,), 0(0),,(0),(0)m m m m m m m d u du d u f t u t dt dt dt d u du v v u v dt dt ----⎧=>⎪⎪⎨⎪===⎪⎩ (i.19) 这时,可以令
1121
(,,,)(,,)m m m du d u u u u u dt dt --=≡u (i.20) 231212(,,,)(,,,,(,,,,))m m m f f f u u u f t u u u ==f (i.21)
012(,,,)m v v v =u (i.22)
于是可以把高阶常微分方程(i.19)化成一阶常微分方程组(i.18)。

i.3收敛性与稳定性
截断误差粗略地说，截断误差可以定义为将微分方程解带入到差分方程后得到的误差，代表了微分方程与差分方程之间的误差。

例如，由Taylor 展式和微分方程(i.1a)得到
221(,())()()()()()(,())2!2!n
n n n n n n n t h h df t u t u t u t hu t u u t hf t u t dt ττ+='''=++=++ 其中n τ是区间1[,]n n t t +上某个常数。

与Euler 法1(,)n n n n u u hf t u +=+相比较，定义余项2(,())2!n
t h df t u t dt τ=为Euler 法的截断误差，它关于h 是2阶的，记为2()h O 。

将上节中讨论的单步法写成一般形式1(,,)0n n F u u h +=，则可以定义截断误差为1((),(),)n n F u t u t h +。

对于每一个差分法在适当的点（例如n t t =，1n t +或者11()2
n n t t ++）作Taylor 展开，就可以得到截断误差的阶。

对多步法可以类似处理。

上节中各方法的截断误差阶分别为：
相容性 一个差分方法称为相容的，如果其截断误差至少是一阶的。

收敛性 实用中我们更关心的是近似解的收敛性，即0h →时，是否有()n n u u t →。

在适当的条件下，例如步长h 足够小，右端函数f 和解u 足够光滑等等，可以证明以上讨论的各方法都是收敛的，并且有估计式
()p n n u t u Ch -≤ （i.23）
其中常数C 与u 和f 有关，与h 和n 无关；阶数p 等于截断误差的阶数减去1。

稳定性 收敛性考虑的是差分方程的精确解与微分方程的精确解之间的误差。

但是在计算机上求差分解时，由于计算机字长的限制，不可避免地产生舍入误差。

另外，方程中的系数和初值等等由于测量条件等限制，也会产生数据误差。

而稳定性研究的就是舍入误差和数据误差对差分方程计算结果的影响，即计算过程中某一步的“差之毫厘”会不会导致后面结果的“失之千里”。

根据误差来源和误差衡量标准的不同，可以定义形形色色的稳定性。

一般说来，隐式方法的稳定性好于同类型的显式方法（例如Euler 方法与隐式Euler 方法）。

另外，为了保证稳定性，常常要求步长h 足够小。

对于常微分方程数值解来说，最简单常用的是关于初值的稳定性。

以单步法为例。

称某差分法关于初值稳定，如果对同一个微分方程和所有足够小的h ，存在常数C ，使得从不同初值0u 和0v 出发的两个差分解{}n u 和{}n v 之间的误差满足
100max n N n n u v C u v ≤≤-≤- （i.24）
其中/N T h =。

由于任意一对n u ，n v 都可以看作是初值，关于初值的稳定性其实是考察如下问题：假设某一步计算有误差，而其后的计算不再有误差，那么这一步的误差对以后结果的影响如何。

本章所讨论的所有差分方法都是关于初值稳定的。

不加区别地说差分方法稳定时，通常指的是关于初值稳定。

不难证明，如果一个差分方法是稳定且相容的，则一定收敛，即
()0, n n u t u -→当n →∞ （i.25）
在差分解的实际计算中，每步都难以避免地舍入误差。

因而，我们进一步考察差分方法的所谓绝对稳定性，即每一步的舍入误差积累起来，是否会对后面的运算结果产生太坏的影响。

考察某一个差分法对非线性常微分方程的绝对稳定性是十分困难的事情，无法得到象对初值稳定性那样的漂亮结论。

通常的做法是考虑模型方程
du u dt
μ= （i.26） 可以认为这个模型方程的真解和差分解的性质代表了非线性常微分方程(,)du f t u dt
=在某一个局部的真解和差分解的性质。

（例如在(,)(,)n n t u t u =处将(,)f t u 线性化，将μ看作(,)n n f t u u
∂∂，忽略其他项。

）记h h μ=。

关于某一个差分法的绝对稳定性的典型结果是：给出一个区间(,)αβ，使得对于所有的(,)h αβ∈，求解（i.25）的这个差分法都是绝对稳定的。

这个区间(,)αβ称之为这个差分法的绝对稳定域。

绝对稳定域越大，h 的取值范围就越大，μ的允许区值就越多，所代表的非线性常微分方程也就越多。

关于绝对稳定性的一个必要条件是：若一个多步法是相容并且稳定的，则它绝对稳定（即它的绝对稳定域非空）的一个必要条件是0μ<。

这相当于仅当
0f u ∂<∂时，多步法才可能绝对稳定。

而当0f u
∂>∂时，在差分计算中误差增长很快。

那么这时候是否就无法用差分法近似计算微分方程了呢？为此，注意到模型方程du u dt
μ=的解是exp()u t μ=。

从而微分方程的解当0μ>（相应于0f u
∂>∂）时是按指数增长的。

因此，如果这时误差的增长速率不超过微分方程解的增长速率，则仍然是可以接受的，并称之为相对稳定。

如果存在区间(,)αβ，使得对于所有的(,)h αβ∈，求解（i.25）的这个差分法都是相对稳定的，则称这个区间(,)αβ为这个差分法的相对稳定域。

特别地，Milne 方法的相对稳定域为(0,)∞。

这就是说，Milne 方法绝对稳定极差，相对稳定性却极好。

习题
写出并解释用于（i.18）的标准Runge-Kutta 方法。