大连理工大学 高等数值分析 常微分方程数值解法-2017

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） 类似地，如果用梯形公式，就得到

111Euler [(,)(,)], 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 法等简单的单步法精度要高一些。下面我们讨论一类所谓Runge-Kutta 法。他们是单步法，但是其精度可以与多步法比美。最常用的是下面的标准Runge-Kutta 法：

标准Lunge-Kutta 法

⎪⎪⎪⎪⎩

⎪⎪⎪⎪⎨⎧++++=++=++=++==+)22(6),()2,2()2,2(),(432113423121K K K K h u u hK u h t f K hK u h t f K hK u h t f K u t f K n n n n n n n n n n （i.12）

从几何上， Runge-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）中的积分。下面介绍Runge-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）