ODE Solver

math

---

我们将分次介绍求解微分方程（ODE）的数值方法。采用的教科书是
Numerical Analysis by Gautschi, Walter, Ch 5 and 6: Initial Value Problems for ODEs.

问题：求解一般的非线性ODE

$$y' = f(t,y),\quad y(a) = y_0.$$
这是典型的初值问题。$t$是一维的自变量，$y = y(t)$是待求的函数，右端 $f(t,y)$是已知的二元函数。假设$f(t,y)$对第二个变量是Lipschitz连续的，即$\|f(t,y_1) - f(t,y_2)\|\leq L\|y_1 - y_2\|$, 则方程的解存在并唯一. 问题是如何数值求解。

除了初值是已知的，从方程我们还可以得知解在初始点的导数，即 $y'(a) = f(a,y_0)$. 由Taylor展示，我们可以得到下一步的逼近

$$y(a+h) \approx y(a) + hy'(a) = y_0 + hf(a,y_0).$$
右边是一个可计算的量，左边是解在 $a+h$ 的值(未知). 为了更好地区分真解和逼近解，我们引入记号 $u_n$ 代表数值解。上述的简单算法可以总结为:
Let $t_0=a,\, u_0 = y_0$, for $n=0,1,2,...$
$$\begin{array} tt_{n+1} &= &t_n + h_n,\\ u_{n+1} &= &u_n + h_n f(t_n,u_n). \end{array}$$
这就是最简单的显示Euler法。

几何上来讲，真解是一条曲线，欧拉法是用折线段来逼近这条解曲线。直观上来看，只要每条边长度逐渐减小，就能更好地逼近真解。数值分析就是精确地分析数值误差，然后改进方法，提高精度。

比如说，欧拉法可以理解成从当前点出发，沿着切线的方向，朝前移动了一小步。有没有更好的方向选择呢？

---

上次介绍了最简单的欧拉法：

$$u_{n+1} = u_n + h_n f(t_n,u_n).$$
欧拉法可以理解成从当前点出发，沿着切线的方向，朝前移动了一小步。有没有比切向更好的前进方向呢？

在继续讨论之前，我们先明确一下好坏的标准。先只考虑一个区间 $[t_n, t_{n+1}]$. 步长记为 $h_n = t_{n+1} - t_n$. 数值解的一般形式是

$$u_{n+1} = u_n + h_n \Phi (t_n, u_n; h_n).$$
其中 $\Phi(t_n, u_n; h_n)$ 可以理解为在数值方法在初始 $(t_n, u_n)$ 处采用的前进方向。显示欧拉法就是 $\Phi(t_n, u_n; h_n) = f(t_n,y_n)$ 的情况。

真解由 ODE $y'=f(t,y)$ 和初值 $y(t_n) = y_n$决定了。真解在 $t_{n+1}$ 的值记为 $y_{n+1} = y(t_{n+1})$. 我们考虑初值相同的简单情况，即 $u_n = y_n$.

好和坏可以由数值解和真解的误差 $u_{n+1} - y_{n+1}$ 来衡量。稍微仔细想一想就知道这个值正比于步长。所以更公平一点的是 平均化的误差

$$T(t_n,y_n;h_n) = \frac{1}{h_n}(u_{n+1} - y_{n+1}).$$
这个称之为 truncation error（截断误差）.

后面的分析可以进一步揭示除上h，实际上是在考虑误差的离散H1范数。Residual 算子在这个范数是适定的.

为了简化记号，我们省略下标 n，下一步的值用下标 next 来替代。这样截断误差可以简写为

$$T(t,y;h) = \frac{1}{h}(u_{\rm next} - y_{\rm next}).$$
我们用下图来总结用到的记号。

根据定义，我们也可以把截断误差写为

$$T(t,y;h) = \Phi (t, y; h) - \frac{1}{h}\int_{t}^{t+h}f(t,y(t)) \rm dt.$$
这个式子让我们可以从另外一个角度来理解截断误差：$\Phi (t, y; h)$ 是对 $f(t,y(t))$ 在区间$[t,t+h]$ 上积分平均的一个近似。

从这个角度出发，我们立马就得到了显示欧拉法的截断误差是一阶，即

$$f(t,y(t))- \frac{1}{h}\int_{t}^{t+h}f(t,y(t)) {\rm d}t = O(h).$$
如果我们用中点公式，或者梯形公式来计算积分，那么误差就是二阶的。更一般的，数值积分可以任意高阶地逼近该积分（高斯点）。有什么额外的困难吗？

我们以中点公式为例。同样由 Taylor 展式，我们容易得到

$$f(t+\frac{h}{2},y(t+\frac{h}{2}))- \frac{1}{h}\int_{t}^{t+h}f(t,y(t)) {\rm d}t = O(h^2).$$
问题在于积分点 $y(t+\frac{h}{2})$ 是未知的，我们只知道初始值 $y(t)$！

怎么办？我们可以先用欧拉法走一半，得到 $u_{\rm half}$, 把它作为 $y(t+\frac{h}{2})$ 的逼近，然后用 $f(t+\frac{h}{2},u_{\rm half})$来替代 $f(t+\frac{h}{2},y(t+\frac{h}{2}))$. 总结起来就是如下的改进的欧拉法

$$u_{\rm next} = u + h f(t+\frac{h}{2},u+\frac{h}{2} f(t,u)).$$
为编程方便，函数 f 的嵌套可以分开如下
$$k_1= f(t,u), \quad k_2 = f(t+\frac{h}{2},u+\frac{h}{2} k_1)$$
$$u_{\rm next} = u + h k_2.$$
换句话说 $k_2$ 是比 $k_1$ 更好的前进方向，截断误差从一阶变成了二阶（请自行证明）。

同样的，对积分的梯形公式做类似的修改，我们得到 ODE solver 的梯形公式，也叫 Heun's 方法.

$$k_1= f(t,u), \quad k_2 = f(t+h,u+hk_1)$$
$$u_{\rm next} = u + h (k_1 + k_2)/2.$$

注意到两种方法都需要算两次斜率。这称之为 two stage 法。更一般的形式是

$$u_{\rm next} = u + (\alpha_1 h) k_1 + (\alpha_2 h) k_2.$$
其中 $\alpha_1\in [0,1),$ 并且 $\alpha_1 + \alpha_2 = 1$。 几何解释是我们沿着 $k_1$ 的方向前进一段距离，然后换个更好的方向，再走完剩下的距离。一般来说 $k_1 = f(t,u)$ 是初始的切线斜率，关键是 $k_2$ 的选取，使得截断误差能得到二阶。