常微分方程（ODE）求解方法总结

1 常微分方程（ODE）介绍

1.1 微分方程介绍和分类

举例：假设跳伞人的下落速度v于时间有如下关系：

$\frac{dv}{dt}=g-\frac{c}{m}v$ (1.1)

其中g为重力常数，m为质量，c为阻力系数。

被微分的量v因变量，与v有关的变量t称为自变量。
如果函数只有一个自变量，那么方程就称为常微分方程（ordinary differential equation，ODE）。
如果函数含有两个或者多个自变量，则成为偏微分方程（partialdifferential equation，PDE）。

此外，微分方程也可以根据阶数来分类：最高阶导数是一阶导数，则称为一阶微分方程（first-order-equation）；最高阶导数是二阶导数，则称为二阶微分方程（second-order-equation）。例如如（1.1）中就是一阶微分方程，下式（1.2）就是一个二阶微分方程。

$m\frac{d^{2}x}{d{t}^2}+c\frac{dx}{dt}+kx=0$ (1.2)

高阶微分方程能简化成一阶方程组。考虑上式（1.2），定义新变量y，令
$y=\frac{dx}{dt}$ (1.3)

对上式取微分得：
$\frac{dy}{dt}=\frac{d^{2}x}{d{t}^2}$ (1.4)

将式（1.3）和（1.4）代入式（1.2）中得到：
$m\frac{dy}{dt}+cy+kx=0$ (1.5)

于是原来的二阶微分方程（1.2）可等价于两个一阶方程组（1.3）和（1.5）。

同样的，其他的n阶微分方程也可以用类似的方式简化。

1.2 常微分方程的非计算机求解方法

常微分方程通常采用解析积分得方法来求解。如对于式(1.1)，先乘以dt，在进行积分得到：
$v=\int (g-\frac{c}{m}v)dt$ (1.6)

对于上式（1.6），是可以精确的推导出该积分得函数表达式的。因为该方程是线性的。

但在实际中，很多方程（是非线性的）精确解是无法求出的。于是提出了一个方法，就是将方程线性化。
（$n$阶）线性常微分方程的一般形式是：
$a_n(x)y^{(n)}+...+a_2(x)y^{(2)}+a_1(x)y^{\prime }++a_0(x)y=f(x)$ (1.7)
其中，$y^{(n)}$是y关于x的n阶导数，$a_n(x)$和$f(x)$都是关于x的函数。因为该方程中未出现因变量y与其导数的乘积，也没有出现非线性函数。所以认为它是线性的。

如下式（1.8）是一个非线性微分方程：
$\frac{d^{2}x}{d{t}^2}+\frac{g}{l}sin(x)=0$ (1.8)
由于含有 $sin(x)$为非线性函数，故该微分方程是非线性的。

线性常微分方程是可以通过解析法求解的。但是，大部分非线性方程无法精确求解。

1.3 线性微分方程求解的推导过程

拿一个简单的方程举例。首先给定函数：
$y=-0.5x^4+4x^3-10x^2+8.5x+1$ （1.9）
这是一个四次多项式。对其进行微分，就得到一个常微分方程：
$\frac{dy}{dx}=-2x^3+12x^2-20x+8.5$ （1.10）

对式（1.10）乘以dx，在进行积分得到：
$y=\int (-2x^3+12x^2-20x+8.5)dx$ （1.11）
应用积分法则得出解为：
$y=-0.5x^4+4x^3-10x^2+8.5x+C$ （1.12）
除了相差一个C外，其余都与原函数相同。这个C称为积分常数（constant of integration）。
出现一个任意常数C表明，积分的结果并不算是唯一的。无限多个常数C对应无限多个可能的函数，都满足微分方程。下图给出了6个满足条件的函数：

为了将解完全确定下来，微分方程通常伴随有辅助条件（auxiliary conditions）。对于一阶常微分方程，有一类被称为**初值（initial value）**的辅助条件，这类条件用于确定常数值，从而使得解是唯一的。例如，给式（1.11）添加初始条件x=0，y=1。带入式（1.12）中，可推导出C=1。于是就得到了唯一解：
$y=-0.5x^4+4x^3-10x^2+8.5x+1$
这个解同时满足常微分方程和指定的初始条件。

当处理n阶微分方程时，就需要n个条件来确定唯一解。如果所有的条件都是在自变量同一值处指定的，那么问题就称为初值问题（initial-value problem）。与之相对的，边值问题（boundary-value problems），就是指在自变量的不同值处指定初始条件。

2 一阶常微分方程（ODE）求解方法

2.1 欧拉方法

2.1.1 欧拉方法

还是拿式（1.1）举例：

$\frac{dv}{dt}\approx \frac{\Delta v}{\Delta t}=\frac{v(t_{i+1})-v(t_i)}{t_{i+1}-t_i}$ （2.1）
其中，$\Delta v$和$\Delta t$分别为速度与时间的差分，$v(t_i)$为初始时刻$t_i$的速度，$v(t_{i+1})$为下一个时刻$t_{i+1}$的速度。注意，$\frac{dv}{dt}\approx \frac{\Delta v}{\Delta t}$ 是一个近似计算，因为：
$\frac{dv}{dt}={\lim }_{\Delta t\to 0}\frac{\Delta v}{\Delta t}$

式(2.1)称为在时刻$t_i$处导数的有限差商(finite divided difference)逼近。将其代入式（1.1）中可得：
$\frac{dv}{dt}\approx \frac{v(t_{i+1})-v(t_i)}{t_{i+1}-t_i}=g-\frac{c}{m}v(t_i)$ （2.2）

对该方程进行整理可得:
$v(t_{i+1})=v(t_i)+[g-\frac{c}{m}v(t_i)](t_{i+1}-t_i)$ （2.3）

注意1：在式（2.3）的中括号内部分$g-\frac{c}{m}v(t_i)$其实就是式（1.1）常微分方程的右边$g-\frac{c}{m}v$，是不过这个是一个估计值，是拿$t_i$对刻的速度来估计$t_i$ ~ $t_{i+1}$这段时间的速度，但只要两个时间点距离越近，这个估值就准确。

根据式（2.3），如果给定了在某一对刻$t_i$的初始值速度$v_i$，就很容易计算出下一时刻$t_{i+1}$的速度$v_{i+1}$。接着又可以用$t_{i+1}$时刻的新速度值来计算$t_{i+2}$的速度，然后依次继续下去。这样，用这种方式就可以计算出任意时刻的速度。
$新值=旧值+斜率\times 步长$
用数学语言表述为:
$v_{i+1}=v_i+\varphi \times h$（2.4）

这种方法的就命名为欧拉法(Euler’s method)(或欧拉一柯西法(Euler-Cauchy method)，或折线法
(point-slope method)。
欧拉方法利用斜率（等于一阶导数在$x_i$处的值)，通过线性外推的办法预测出y在前进步长h之后$x_{i+1}$处的新值：

利用上面的方法，取步长为2（也就是每隔2秒，计算一次速度值），于是得到一系列时间t与速度v的值，将这些计算结果与精确解同时画在一个图中：

图5 速度v与时间t的关系的常微分方程的数值解与解析解的对比

上图中，如果我们将步长减小一半，以1秒的步长计算得到的结果，其误差将会更小，从而直线段的轨迹更趋近于真实解。但是步长每减小一半，可以得到更精确的结果，所需要的计算量也会翻一番。

2.1.2 欧拉方法的误差分析

常微分方程数值解的误差包括两种类型：
(1)截断(Truncation)误差或离散误差，是由逼近Y值的方法本身引起的。
(2)舍入(Round-off误差，是由于计算机能保留的有效数字有限而引起的。
截断误差由两部分组成：第一部分是局部截断误差(local truncation error)，指方法在一步迭代中所产生的误差。第二部分是传播截断误差(propagated truncation error)，它是由前面各步中产生的逼近值引起的，这两部分的和称为整体截断误差或全局截断误差(global truncation error)。

利用泰勒级数展开直接推导欧拉方法，可以同时推出截断误差的大小和属性。假设被积分的微分方程具有一般形式：
$y^{\prime }=f(x,y)$（2.5）
其中，$y^{\prime }=dy/dx$，x和y分别为自变量和因变量。

如果解(也就是描述y行为的函数)具有连续导数，那么它可以在点(xi, yi)处展开成泰勒级数，即
$y_{i+1}=y_i+y_i^{\prime }\ast h+\frac{y_i^{\prime \prime }}{2!}\ast h^2+...+\frac{y_i^{(n)}}{n!}\ast h^n+R_n$（2.6）
其中，$h=(x_{i+1}-x_i)$，$R_n$为余项，定义为：
$R_n=\frac{y^{(n+1)}(\xi )}{(n+1)!}\ast (h{)}^{(n+1)}$（2.7）
其中，$\xi \in [x_i,x_{i+1}]$。

由于$y^{\prime }=f(x,y)$，所以$y^{\prime \prime }=f^{\prime }(x,y)$，$y^{\prime \prime \prime }=f^{\prime \prime }(x,y)$，依此类推，代入式（2.6）右边，同时也加入式（2.7），得到：
$y_{i+1}=y_i+f(x_i,y_i)\ast h+\frac{f^{\prime }(x_i,y_i)}{2!}\ast h^2+...+\frac{f^{(n-1)}(x_i,y_i)}{n!}\ast h^n+O(h^{n+1})$（2.8）
此处$O(h^{n+1})$表示局部截断误差与步长的n+1次幂成正比。

因为我们仅使用了泰勒级数的有限项来逼近真解，所以就产生了截断误差。这样，我们就截断或者略去了一部分真实解。例如，在前面的计算中，我们采用$y_{i+1}=y_i+f(x_i,y_i)\times h$来计算下一个点的值，有部分泰勒级数展开项没有被包括在式中，这些剩余项就导致了欧拉方法的截断误差：
$E_t=\frac{f^{\prime }(x_i,y_i)}{2!}\ast h^2+...+\frac{f^{(n-1)}(x_i,y_i)}{n!}\ast (h{)}^n+O(h^{n+1})$（2.9）
$E_t$为真实的局部截断误差。

当h足够小时，式(2.9)中各项的误差通常随着阶数的增加而减少。于是常常可以表示为：
$E_a=\frac{f^{\prime }(x_i,y_i)}{2!}\ast h^2$（2.10）
或者：
$E_a=O(h^2)$（2.11）
$E_t$为近似的局部截断误差。

泰勒级数给出了一种量化欧拉方程误差的方法。不过，这种方法使用时存在局限:
(1)泰勒级数仅能估计局部截断误差一即方法在一步迭代中所产生的误差。它不能估量传播截断误差，因此也不能估量全局截断误差。
(2)实际问题中处理的函数一般比简单多项式复杂。因此，计算泰勃级数展开时，导数值有时很难得到。

尽管这些限制使泰勒级数无法分析人多数实际问题的精确误差，但是它仍然有助于了解欧拉方法的性态。按照式(2.11)可发现局部误差与步长的平方和微分方程的一阶导数成正比。此外已经证明，全局截断误差为O(h)，全局截断误差与步长成正比(Carnahan等人，1969)。

结论:
(1)缩短步长可以减少误差。
(2)因为直线的二阶导数为常数，所以如果基本函数（即微分方程的解)是线性的，那么方法能给出无误差的预测值。

由于欧拉方法使用直线段来逼近解，因此，欧拉方法称为一阶方法(first-order method)。
思考：如果使用二次来逼近解，也就是二阶方法，就能给出基本函数（即微分方程的解)是二次的问题的精确。同理，如果基本函数是n次多项式，则使用n阶方法可得出精确解。

2.1.3 欧拉方法的改进思路1——添加高阶项

减小欧拉方法的误差的一种方法是，在公式中加入泰勒级数展开的高阶项。例如，添加二阶项后，得到：
$y_{i+1}=y_i+f(x_i,y_i)\ast h+\frac{f^{\prime }(x_i,y_i)}{2!}\ast h^2$（2.12）
于是局部截断误差变为：
$E_a=\frac{f^{\prime \prime }(x_i,y_i)}{3!}\ast h^3$（2.13）

尽管添加高阶项对于处理多项式来说十分简单，不过有时候会使常微分方程变得更加复杂。特别地，如果常微分方程同时为自变量和因变量的函数，就需要通过链式法则来计算微分。例如，f(x,y)的一阶导数是：
$f^{\prime }(x,y)=\frac{\partial f(x_i,y_i)}{\partial x}+\frac{\partial f(x_i,y_i)}{\partial y}\frac{dy}{dx}$
二阶导数是：
$f^{\prime \prime }(x,y)=\frac{\partial [\partial f/\partial x+(\partial f/\partial y)(dy/dx)]}{\partial x}+\frac{\partial [\partial f/\partial x+(\partial f/\partial y)(dy/dx)]}{\partial y}\frac{dy}{dx}$

高阶导数会变得越来越复杂。而综合其计算量，其在性能上其实与一阶泰勒级数方法差不多，而一阶泰勒展开的欧拉方法只需要计算一阶导数，实现起来会简洁很多。因此，添加高阶项的方法其实并不是很实用。

2.1.4 欧拉方法的改进思路2——对斜率采用更好的估值方法

在“注意1”中我们提到，欧拉方法是将区间左端点处的斜率当作斜率在整个区间上平均值的逼近。但估值毕竟是估值，不是精确值，所以这个估值的精确度会严重影响计算的结果。
于是我们思考，有没有什么其他估计斜率的其他方法，能够给出更好的估值，从而得到更精确的预测值。
基于这个想法，下面给出了两种简单的修正方法来克服这个估值不精确的问题：修恩法和中点方法。
(这两种方法其实都属于一个更广泛的类型—龙格一库塔法。)

2.1.4.1 修恩法

改进斜率估计值的一种方法是，在区间上计算两个导数值：一个在区间的起始点处，一个在区间的终点处。然后，将这两个导数的平均值作为斜率在整个区间上的改进估计值。这种方法，称为修恩法(Heun’s method)，其示意图如下图所示。

图5 修恩法的示意图：(a)预估和(b)校正

前面介绍过，在欧拉方法中，区间起始点处的斜率为：
$y_i^{\prime }=f(x_i,y_i)$ (2.15)
利用这个值线性外推出$y_{i+1}$：
$y_{i+1}^0=y_i+f(x_i,y_i)h$ (2.16)
标准的欧拉方法算到这一步就停止了。然而，在修恩法中，由式(2.16)得出的$y_{i+1}^0$并不是最终结果，而是一个中间的预测值。式(2.16)称为预估表达式(predictor equation)。该式给出了$y_{i+1}$的一个估计值，利用这个估计值又可以计算区间终点处的斜率（形式如式(2.15)）：
$y_{i+1}^{\prime }=f(x_{i+1},y_{i+1}^0)$ (2.17)
我们得到两个斜率估计值$y_i^{\prime }$和$y_{i+1}^{\prime }$(式(2.15)和式(2.17))，于是可以计算斜率的区间平均值：
${\overline{y_i}}^{\prime }=\frac{y_i^{\prime }+y_{i+1}^{\prime }}{2}=\frac{f(x_i,y_i)+f(x_{i+1},y_{i+1}^0)}{2}$ (2.18)

然后，按照欧拉方法，利用这个平均斜率由$y_i$线性外推出$y_{i+1}$:
$y_{i+1}=y_i+\frac{f(x_i,y_i)+f(x_{i+1},y_{i+1}^0)}{2}h$ (2.19)
此式称为校正表达式((corrector equation)。
修恩法是一种预估一校正方法{predictor-corrector approach）。

修恩法可简洁地表示为：
预估：$y_{i+1}^0=y_i+f(x_i,y_i)h$ (2.16)
校正：$y_{i+1}=y_i+\frac{f(x_i,y_i)+f(x_{i+1},y_{i+1}^0)}{2}h$ (2.19)

迭代修恩法：
因为式(2.19)的等号两边都含有$y_{i+1}$，所以该式可以按迭代方式计算，可以反复使用旧的估计值计算出改进后的$y_{i+1}$值。// 需要注意的是，这个迭代过程并不一定能收敛到真实解，不过它将收敛到一个截断误差有限的估计值。

校正过程收敛的终止条件为：
$|{\epsilon }_a|=|\frac{y_{i+1}^j-y_{i+1}^{j-1}}{y_{i+1}^j}|100\%$
其中，$y_{i+1}^{j-1}$和$y_{i+1}^j$分别为上一校正步的结果和当前校正步的结果。

// 前面没有迭代校正的修恩法，可称为简单修恩法；经过迭代计算校正的修恩法，可以称为迭代修恩法。

特殊情况：
在前面的例子中，导数$dy/dx$同时为因变量y和自变量x的函数（$dy/dx=f(x,y)$）。对于常微分方程仅与自变量有关的情况（$dy/dx=f(x)$），求解过程中并不需要预估步（式(2.16)），而且每次迭代中也只需一次校正。在这种情况下，方法简化为
$y_{i+1}=y_i+\frac{f(x_i)+f(x_{i+1})}{2}h$ (2.19)

该式与“数值积分方法的总结”第1.1.1节 梯形法则（2点积分）其实是等价的。相应的截断误差也与梯形法则一样为：
$E_t=-\frac{f^{\prime \prime }(\xi )}{12}{h}^3$

由于真解是二次多项式时，常微分方程的二阶导数为0，所以方法是二阶精度的。另外，局部和全局误差分别为$O(h^3)$和$O(h^2)$，因此，减小步长时，误差减小的速度比欧拉方法快。

2.1.4.2 中点方法

图5和式(2.16)（$y_{i+1}^0=y_i+f(x_i,y_i)h$ ）使用了端点处的值$y_{i+1}$ 来修正斜率的。

那么是不是也可以考虑用中点处的值来修正斜率？？

于是产生了另一种简单的欧拉方法的修正方法（如图6所示）。该方法称为中点方法(midpoint method)(或改进的多边形法（improved polygon），或修正的欧拉方法(modified Euler))，利用欧拉方法预测区间中点处的y值（图6(a)）：
$y_{i+1/2}=y_i+f(x_i,y_i)h/2$ (2.20)
然后用这个值来计算中点处的斜率：
$y_{i+1/2}^{\prime }=f(x_{i+1/2},y_{i+1/2})$ (2.21)
假设这个值是对整个区间的斜率平均值的有效近似，那么由$x_i$线性外推出$x_{i+1}$处的值为:
$y_{i+1}=y_i+f(x_{i+1/2},y_{i+1/2})h$ (2.22)

// 注意对比式(2.20)、(2.20)、(2.20)与 (2.16)、 (2.17)、 (2.19)，就可以清晰地看出中点法和修恩法的区别。

图6 中点公式的示意图。(a)式(2.20)}和(b)式(2.22)

因为式（2.22）只有等式的一边含$y_{i+1}$，所以这种方法的结果无法像修恩法一样采用迭代校正的办法来改进。

（本方法也可同牛顿一柯特斯开型公式的中点公式结合起来理解。）

因为中点方法利用区间中点处的斜率估计值，它理论上比端点（起点）处的斜率估计值要更好(//关于数值微分的讨论知，中心有限差分逼近导数的效果比前向或后向差分都要好)，所以中点方法优于欧拉方法。中点方法的局部和全局误差分别为$O(h^3)$和$O(h^2)$。

2.1.5 欧拉方法的改进思路3

思考：上一节提到的中点法，就类似计算积分的时候，取中点计算估值（1点积分）；修恩法，就类似于梯形法计算积分（2点积分）。回想在“数值积分方法的总结”积分计算中，除了梯形法则（2点积分），我们还讲到了3点积分和4点积分，区别就在于用二次曲线和三次曲线来逼近使积分结果更加接近真实值。那么，这个地方，如果我们能计算或者估算出三个点或者四个点，是否也能用类似的方法，使求解的结果更接近真实值？

2.2 龙格-库塔法

龙格一库塔法(Range-Kutta)(RK)无需计算高阶导数，也可以达到与泰勒级数方法同样的精度。龙格一库塔法存在许多不同版本，但它们都可以写成式(25.1)的一般形式
$y_{i+1}=y_i+\varphi (x_i,y_i,h)h$ (2.23)
其中，$\varphi (x_i,y_i,h)$称为增量函数((increment function)，它表示整个区间上的斜率。增量函数的一般形式是

k1出现在k2的表达式中，k2又出现在k3的表达式中，等等。因为每个k都是由函数赋值的，所以这种
递归关系提高了龙格-库塔法的计算效率。

将增量函数的项数记为n，n值的不同就对应着不同类型的龙格一库塔法。n=1时，一阶龙格-库塔法，实际上就是欧拉方法。

（未完待续，更新中。。。）