大连理工大学 高等数值分析 常微分方程数值解法-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 ），得11()()(,())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(,)()2n 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 ），得22()()(,())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）选取这些待定常数,,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 可以依次显式地算出。

如果在某一个i K 的表达式中出现j K ，其中j i >，则导致隐式Runge-Kutta 方法，可以迭代求解。

一般来说，隐式Runge-Kutta 方法的稳定性更好一些。

i.2 常微分方程组与高阶常微分方程先来考虑下面的常微分方程组初值问题11122212121100(,,,,)(,,,,)(,,,,)(0),,(0)m m mm m m m du f t u u u dt du f t u u u dt du f t u u u dtu 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!nn 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!nt 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()2n n t t ++）作Taylor 展开，就可以得到截断误差的阶。

对多步法可以类似处理。

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

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

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