【高等数值分析】常微分方程数值解

1. 预备理论

求解常微分方程初值问题数值解 \[ \begin{align} &\frac{dy}{dx} = f(x,y), \quad a < x < b, |y| < \infty \\ &y(a) = y_0 \end{align} \] 存在唯一性定理：若 \(f(x,y)\) 连续，对 \(y\) 满足 Lipschitz 条件，那么初值问题有唯一解。

对上面的常微分方程积分就有 \(y(t_{n+k}) - y(t_{n-j}) = \int_{t_{n-j}}^{t_{n+k}} f(t,y(t)) dt\)，所以实际上可以转化为数值积分的问题，所以这一节的方法和数值积分很类似。但是这里的问题在于被积函数值是不知道的。

2. 线性多步法

在数值积分里面 \(f(t_k,y(t_k))\) 是已知的，主要是在设计求积节点对应的系数。在这里还需要首先估计每个节点的函数值。对前面的式子采样求积就得到 \(y(t_{n+k}) - y(t_{n-j}) = h\sum_i \beta_i f_{n-i}\)，就引出了下面要介绍的线性多步法。

2.1 线性多步法

记 \(f_n=f(t_n,y_n), y_n=y(t_n)\)，线性多步法（线性 \(k\) 步法）的一般表达式为 \[ \sum_{i=0}^k \alpha_i y_{n+i} = h\sum_{i=0}^{k} \beta_i f_{n+i} \] 其中 \(\alpha_k=1, \alpha_0^2+\beta_0^2\ne0\)。给一些特例：

• \(\alpha_0=-1,\beta_0=1,\beta_1=0\) 时为显式Euler法 \(y_{n+1}=y_n + hf_n\)；
• \(\alpha_0=-1,\beta_0=0,\beta_1=1\) 时为隐式Euler法 \(y_{n+1}=y_n + hf_{n+1}\)；
• \(\alpha_0=-1,\beta_0=\beta_1=1/2\) 时为梯形法 \(y_{n+1} = y_n + h/2(f_n+f_{n+1})\)。

可以定义“特征多项式”，在后面分析稳定性和收敛性的时候会很有用 \[ \rho(\xi) = \sum_{i=0}^k \alpha_i \xi^i, \quad \sigma(\xi)=\sum_{i=0}^k \beta_i\xi^i \] 如何衡量数值求解方法的精度呢？定义局部截断误差为 \[ T_{n+k} = \sum_{i=0}^k [\alpha_i y(t_{n+i}) - h\beta_i y'(t_{n+i})] \] 对上面的式子做Taylor展开，可以得到形如 \(T_{n+1} = -1/90 h^5 y^{(5)}(x_n)+O(h^6)\)（这是Simpson公式的局部截断误差），这个式子表明Simpson公式是4阶的。

对于线性 \(k\) 步法，要想设计一种迭代方法使得数值精度是 \(p\) 阶的，可以采用待定系数法，保证局部截断误差中 \(y^{(1)}(x),...,y^{(q)}(x)\) 前面的系数都是0。

2.2 稳定性与收敛性

2.2.1 相容性

首先引入相容性的概念。若某个多步法的截断误差满足 \(T(t,y(t),h)=o(h)\) 那么称其为相容的。若 \(T(t,y(t),h)=O(h^{q+1})\)，则称其为 \(q\) 阶的。

实际上，相容性就是说该方法至少是1阶的，也就是当 \(y\) 为一次线性函数时，该方法要能够得到准确解。

定理：相容的多步法充要条件是 \(\rho(1)=0,\rho'(1)=\sigma(1)\)。

2.2.2 零稳定性

根条件：若多项式 \(\rho(\xi)\) 的 \(k\) 个根的模长都不大于 1，并且模值等于 1 的根都是单根，则称其满足根条件。

Note：根条件考虑的是齐次差分方程的解。当 \(f\equiv 0\) 的时候，齐次差分方程的解具有指数形式 \(y(n)=\sum_i P(n) \xi_i^n\)，其中 \(P(n)\) 是一个多项式， \(\xi_i\) 就是 \(\rho(\xi)\) 的根（这里没考虑重根的情况，不过是类似的），如果 \(|\xi_i|<1\) 就意味着 \(\lim_{n\to\infty}y(n) = 0\)，那么这跟稳定性有什么关系呢？

定理：线性多步法关于初值稳定的充要条件是 \(\rho(\xi)\) 满足根条件。

Note：这个定理说明了根条件和（零）稳定性二者的等价关系。为什么呢？假设真实的初值为 \(y_0=y(t_0)\)，真实的解为 \(y(n)=y_0\xi_0^{(n-t_0)}\)。而我们计算的时候由于各种原因拿到的初值是有误差的，也就是 \(\hat{y}_0\)，那么最后求得的误差就是 \(\hat{y}(n)-y(n)=(\hat{y}_0-y_0)\xi_0^{(n-t_0)}\)，如果不满足根条件，当 \(n\to\infty\) 的时候，\(\hat{y}(n)-y(n) \to \infty\)，解就是不稳定的。

2.2.3 收敛性

定义：假设在区间 \([a,b]\) 上等距划分 \(N\) 个区间，\(x_n = a+nh,n=0,1,...,N\)，求解初值问题得到的解为 \(y_n,n=0,1,...,N\)，最大整体误差为 \(E(h) = \max_{n} |y(x_n) - y_n|\)，如果满足 \(\lim_{h\to0} E(h)=0\)，则称方法是收敛的。如果 \(E(h) \le Ch^p\)，则称方法是 \(p\) 阶收敛的。

定理：相容性 + 零稳定性 \(\iff\) 收敛性。

2.2.4 绝对稳定性

前面的零稳定性是跟所选择的方法有关的。而这里的绝对稳定性研究的是微分方程本身的属性。如果方程本身性质很不好，那么可能无论选择什么方法都是不稳定的。

取 \(f(x,y)=\lambda y\)，再研究方法的稳定性。以试验方程为例 \[ \begin{cases} y' = \lambda y, t\in[a,b] \\ y(a) = \tilde{y}_0 \end{cases} \] 用Euler法得到 \(y_n = (1+h\lambda)y_{n-1}\)，于是两个解的误差满足 \(y_n-z_n = (1+h\lambda)^n (y_0-z_0)\)，若 \(|1+h\lambda|>1\)，则误差总会放大，这是我们不希望的。

要保证稳定性，既与方程本身的性质（也即 \(\lambda\)） 有关，也与所选择的步长 \(h\) 有关。从另一个角度而言，方程本身的性质（\(\lambda\)）会影响可选择的 \(h\) 的范围。由此引出绝对稳定性的概念。

对前面提到的线性多步法，把 \(f(x,y)=\lambda y\) 代回去就有 \[ \sum_{i=0}^k (\alpha_i - h\lambda \beta_i) y_{n+i} = 0 \] 称 \(\Pi(\xi;z) = \rho(\xi) - z\sigma(\xi)\) 为稳定多项式。

定理：对给定的 \(z\)，若 \(\Pi(\xi;z)\) 的 \(k\) 个根的模都小于 1，则其是绝对稳定的。（满足这样条件的 \(z\) 的集合构成绝对稳定区域。）

3. Runge-Kutta方法

Runge-Kutta方法和线性多步法的主要区别在于，其在 \(x_n,x_{n+1}\) 中间又进行了采样、插值。这种方法大概可以表示为 \[ \begin{align} & y_{n+1} = y_n + h\sum_{i=1}^{s} b_i k_i \\ & k_i = y_n + f(t_n + c_i h, y_n + \sum_{i=1}^s a_{ij} k_j) \end{align} \] 其中 \(0\le c_i\le1\)。实际上 \(k_i \approx y(t_n+c_i h)\)，就是在 \([y_n,y_{n+1}]\) 中间又进行了采样插值。

这样提高了精度，同时也会增加计算复杂度。

除此之外，提高精度的方法还有 Richardson 外推方法。

4. 刚性问题

刚性问题也是方程本身的属性，主要是指某些情况下两个特解的尺度相差很大，比如两个指数衰减的过程混合在一起，但是其中一个衰减特别快（\(\exp(-\lambda t)\) 的 \(\lambda\) 特别大），另一个则衰减特别慢，那么数值求解的时候很可能只能看到衰减慢的那个过程，另一个则被忽略。这在化学反应中是经常遇到的。