DLSODA：自动切换 Adams 与 BDF 的常微分方程求解器（一）

[TOC]

DLSODA 是 ODEPACK 中最著名的常微分方程初值问题求解器之一。如果要用一句话概括它的特点，那就是：

从非刚性的隐式 Adams 方法开始，在积分过程中判断 Adams 与 BDF 哪一种更有效，并自动完成方法、阶数和步长的切换。

这里的 D 表示 double precision，LSOD 来自 Livermore Solver for Ordinary Differential Equations，最后的 A 表示 automatic method selection。通常文献中讨论算法时写作 LSODA，调用双精度 Fortran 子程序时则写作 DLSODA。

DLSODA 实际上是一个及其复杂的求解器，它包含以下组成部分：

1. 线性多步离散格式；
2. Nordsieck 历史数组和预测&校正过程；
3. 局部误差估计、变步长和变阶；
4. 函数迭代与 chord/Newton 迭代；
5. Jacobian 的差分近似、复用和 LU 分解；
6. 根据“另一种方法可能采用的步长”进行自动方法切换；
7. 对舍入误差、失败重试和输出插值的工程化处理。

本文从这些部分出发，解释 DLSODA 背后的数学原理，并展示其程序实现。

---

1. DLSODA 求解什么问题

DLSODA 处理显式一阶常微分方程组的初值问题

$$\frac{d\mathbf y}{dt}=\mathbf f(t,\mathbf y), \qquad \mathbf y(t_0)=\mathbf y_0,$$

其中

$$\mathbf y(t)\in \mathbb R^N, \qquad \mathbf f:\mathbb R\times\mathbb R^N\rightarrow\mathbb R^N.$$

需要注意的是，它不是一般的微分代数方程（DAE）求解器。若问题写成

$$A(t,\mathbf y)\dot{\mathbf y}=\mathbf g(t,\mathbf y),$$

且矩阵 可能奇异，则应考虑 ODEPACK 中的 LSODI/LSODIS，或者现代 DAE 求解器，而不能直接把问题交给 DLSODA。

DLSODA 的两个积分分支是：

问题的局部性质	积分方法	默认最高阶	非线性迭代
非刚性	隐式 Adams 方法	12	函数迭代
刚性	BDF 方法	5	chord/Newton 迭代

这里我们提到所谓问题的“局部性质“，这是因为同一个问题可以在某些时间区间表现为非刚性，在另一些区间表现为刚性。DLSODA 允许在一次积分中来回切换，这也是这个求解器最大的特色。

---

2. 刚性方程（stiff ODE）

为了理解什么是刚性方程，我们先来看看两个简单的例子

2.1 两个简单的例子

考察最简单的衰减方程

$$y'=-y,\qquad y(0)=1,$$

精确解为

$$y(t)=e^{-t}.$$

使用显式欧拉法：

$$y_{n+1}=(1-h)y_n.$$

从保持数值解的稳定性（与收敛性的概念略有不同）方面考虑，我们对 有如下要求：

$$|1-h|<1,$$

即

$$0<h<2.$$

但是这并不代表只要选择这个区间内的 $h$ 就可以得到精确的结果。从精度的角度考虑，对于这样一个简单的例子，我们可以计算出显式欧拉法在终点处的误差为

$$E(h)=e^{-T}-(1-h)^{T/h}.$$

经过简单的泰勒展开分析之后，我们有结论

$$E(h)=O(h).$$

也就是说，越小的步长对应着越小的数值误差，这跟我们的常识是一致的。

对于 的情形，很显然，在选择步长的时候精度的限制要强于稳定性的限制。这种由精度决定步长的方程，我们称之为非刚性方程。

下面我们再考虑一个更加复杂一点的方程，

$$y'=-1000(y-\cos t)-\sin t,\qquad y(0)=1.$$

这个方程的精确解是

$$y(t)=\cos t.$$

显式欧拉法从精确值出发，一步产生的局部误差为

$$d_{n+1} = \cos t_n-h\sin t_n-\cos(t_n+h) = \frac{h^2}{2}\cos t_n+O(h^3).$$

假设局部误差容限为

$$\tau=10^{-4},$$

那么精度条件近似为

$$\frac{h^2}{2}\le10^{-4},$$

即

$$h\le\sqrt{2\times10^{-4}} \approx0.014.$$

因此

$$h_{\mathrm{acc}}\approx0.014.$$

现在我们再来考察其稳定性的要求：

$$|1-1000h|\le1,$$

也就是说

$$h_{\mathrm{stab}}=0.002.$$

可以看到，稳定性对于步长的要求要大于精度的要求。这种稳定性要求决定步长的方程，我们就称之为刚性方程

总的来说，在ODE的数值计算中，“刚性方程”通常指

方程的解中同时存在差异很大的时间尺度，使得显式数值方法为了保持稳定，不得不采用远小于实际精度需求的步长。

2.2 刚性比

在线性化意义下，令

$$J(t,\mathbf y)=\frac{\partial \mathbf f}{\partial \mathbf y}$$

为 Jacobian。若 的特征值具有相差很大的负实部，例如

$$|\operatorname{Re}\lambda_{\max}| \gg |\operatorname{Re}\lambda_{\min}|,$$

系统通常同时包含快速衰减与缓慢演化的模态。由此我们可以定义刚性比

$$\kappa_{\mathrm{stiff}} = \frac{\max_i|\operatorname{Re}\lambda_i|} {\min_i|\operatorname{Re}\lambda_i|}$$

来表征方程“刚化”的程度。不过它的定义并不严格：Jacobian 随时间变化、非正规矩阵、振荡模态和非线性效应都会使单一特征值比值失效。

当比较大的时候，方程的解就会同时存在快速衰减和缓慢变化的分量，而步长受限于快速衰减分量的要求。这就使得当快速衰减分量衰减完毕之后，我们仍然使用很小的步长来求解缓慢变化的分量，导致求解效率下降。

2.3 LSODA的解决方案

LSODA对于刚性方程的处理方案非常直接，简单来说就是

自动判断当前问题是否表现出刚性，并在非刚性方法与刚性方法之间自动切换。

具体地，LSODA使用了以下两类方法（两种方法的具体内容将在下一章节介绍）：

问题状态	使用的方法	特点
非刚性	Adams	显式预测—隐式校正，单步成本较低
刚性	BDF	隐式方法，稳定区域大，适合刚性问题

LSODA 会监测求解过程中的信息，例如：

• 当前方法能够稳定采用的步长；
• 误差估计和步长调整情况；
• 校正迭代是否收敛；
• 使用另一类方法预计能否取得更高效率。

如果它发现 Adams 方法需要使用远小于精度要求的步长，并判断 BDF 方法会更有效，就会认为问题出现了刚性，并切换到 BDF。反过来，如果刚性消失、Adams 方法预计更高效，LSODA也可以从 BDF 切回 Adams。

下面我们将详细介绍LSODA具体是如何在两种方法之间来回切换的。

---

3. 统一框架：线性多步法

不论是Adams方法还是BDF方法，它们都可以用统一的框架来描述——k步线性多步法。

3.1 一般形式

步线性多步法可以写成

$$\sum_{j=0}^{k}\alpha_j\mathbf y_{n+j} = h\sum_{j=0}^{k}\beta_j \mathbf f_{n+j},$$

其中

$$\mathbf f_{n+j} =\mathbf f(t_{n+j},\mathbf y_{n+j}).$$

如果 ，右端包含未知的 ，方法就是隐式的，需要在每一步求解非线性方程。

这里的系数 和 并不是任意选取的，具体要求可以参考讨论k步线性多步法的相关文章，这里不再赘述。但是不同的选择就对应着不同的数值方法，它们的求解精度也不尽相同。

DLSODA 选用两类隐式线性多步法：

• Adams-Moulton 型方法：高阶精度好，适合非刚性问题；
• Backward Differentiation Formula（BDF）：稳定性更强，适合刚性问题。

3.2 Adams 方法

Adams 方法的出发点是对 插值，然后积分

$$\mathbf y_{n+1} = \mathbf y_n+ \int_{t_n}^{t_{n+1}}\mathbf f(t,\mathbf y(t))\,dt.$$

引入无量纲坐标

$$t=t_n+\theta h,\qquad 0\leq\theta\leq1,$$

则

$$y(t_{n+1}) = y(t_n)+h\int_0^1 f(t_n+\theta h,y(t_n+\theta h))\,d\theta.$$

因此问题变成：

$$\text{如何用已知的离散函数值近似 } \int_0^1 f(t_n+\theta h)\,d\theta?$$

它的一种方法是Adams-Bashforth方法，它是一种显式方法。对于s步的Adams-Bashforth方法，我们使用s个历史斜率

$$f_n,f_{n-1},\ldots,f_{n-s+1}.$$

它们对应于无量纲节点

$$\theta=0,-1,\ldots,-(s-1).$$

用这些节点构造 (s-1) 次Lagrange插值多项式

$$P_{s-1}(\theta) = \sum_{j=0}^{s-1}f_{n-j}L_j(\theta),$$

其中

$$L_j(\theta) = \prod_{\substack{m=0\\m\neq j}}^{s-1} \frac{\theta+m}{m-j}.$$

将插值多项式积分：

$$y_{n+1} = y_n+h\int_0^1P_{s-1}(\theta)\,d\theta.$$

于是得到任意 (s) 步 Adams–Bashforth 公式

$$\boxed{y_{n+1} = y_n+h\sum_{j=0}^{s-1}b_j^{(s)}f_{n-j}}$$

其中系数的一般表达式是

$$b_j^{(s)} = \int_0^1 \prod_{\substack{m=0\\m\neq j}}^{s-1} \frac{\theta+m}{m-j}\,d\theta$$

这些系数完全由插值节点决定。

前几阶为：

$$\begin{aligned} s=1:\quad y_{n+1} &=y_n+hf_n, \\[2mm] s=2:\quad y_{n+1} &=y_n+h\left( \frac32f_n-\frac12f_{n-1} \right), \\[2mm] s=3:\quad y_{n+1} &=y_n+h\left( \frac{23}{12}f_n -\frac{16}{12}f_{n-1} +\frac{5}{12}f_{n-2} \right), \\[2mm] s=4:\quad y_{n+1} &=y_n+h\left( \frac{55}{24}f_n -\frac{59}{24}f_{n-1} +\frac{37}{24}f_{n-2} -\frac{9}{24}f_{n-3} \right). \end{aligned}$$

它的另一种方法是Adams-Moulton方法，它是一种隐式方法。类似的，我们使用s个历史斜率

$$f_{n+1},f_n,f_{n-1},\ldots,f_{n-s+1},$$

对应节点

$$\theta=1,0,-1,\ldots,-(s-1).$$

为了统一编号，令

$$r=-1,0,\ldots,s-1,$$

我们利用Lagrange差值多项式来近似计算积分，

$$P_{s}(\theta) = \sum_{r=-1}^{s-1}f_{n-r}\widetilde L_r(\theta),$$

其中

$$\widetilde L_r(\theta) = \prod_{\substack{m=-1\\m\neq r}}^{s-1} \frac{\theta+m}{m-r}.$$

这样我们就可以把积分写为

$$y_{n+1} = y_n+h\int_0^1P_{s}(\theta)\,d\theta.$$

将多项式的具体表达式带入，我们就可以得到s步的Adams-Moulton公式

$$\boxed{y_{n+1} = y_n+ h\sum_{r=-1}^{s-1}a_r^{(s)}f_{n-r}}$$

其中

$$a_r^{(s)} = \int_0^1 \prod_{\substack{m=-1\\m\neq r}}^{s-1} \frac{\theta+m}{m-r}\,d\theta$$

特别地，(a_{-1}^{(s)}) 是 (f_{n+1}) 的系数。

比方说对于二阶隐式Adams，其刚好就是梯形公式

$$y_{n+1} = y_n+\frac h2(f_{n+1}+f_n)$$

对于三阶Adams-Moulton为

$$y_{n+1} = y_n+h\left( \frac5{12}f_{n+1} +\frac23f_n -\frac1{12}f_{n-1} \right)$$

正常来说，隐式方法的求解涉及到非线性方程的求解，我们一般会采用迭代法。对于牛顿迭代法而言，我们需要计算和分解方程的雅可比矩阵。但对于非刚性问题，这通常并不划算，因此LSODA的Adams分支倾向采用便宜的校正迭代。

3.3 Adams方法的预测-校正形式

在LSOD中，我们经常把两种Adams方法配合起来：

1. 用 (s) 阶 Adams–Bashforth 预测：

$$y_{n+1}^{(0)} = y_n+h\sum_{j=0}^{s-1}b_jf_{n-j};$$

1. 计算预测斜率：

$$f_{n+1}^{(0)} = f(t_{n+1},y_{n+1}^{(0)});$$

1. 用 (s+1) 阶 Adams–Moulton 校正：

$$y_{n+1}^{(1)} = y_n+h\left[ a_{-1}f_{n+1}^{(0)} +\sum_{j=0}^{s-1}a_jf_{n-j} \right];$$

1. 必要时继续迭代。

这就是常见的 PECE 结构：

$$\boxed{ \text{Predict} \rightarrow \text{Evaluate} \rightarrow \text{Correct} \rightarrow \text{Evaluate} }$$

预测值还可以提供很好的隐式迭代初值。

3.4 BDF 方法

BDF 的思路是用当前未知解 (y_{n+1}) 和若干历史解构造插值多项式来预测在最新时刻 (t_{n+1}) 的导数，并令该导数等于 (f(t_{n+1},y_{n+1}))。

具体地，考虑

$$y'(t)=f(t,y(t)), \qquad t_n=t_0+nh.$$

为了从 (t_n) 推进到 (t_{n+1})，取 (k+1) 个解值

$$y_{n+1},y_n,y_{n-1},\ldots,y_{n+1-k}.$$

构造一个 (k) 次插值多项式 (P_k(t))，使

$$P_k(t_{n+1-j})=y_{n+1-j}, \qquad j=0,\ldots,k.$$

然后用

$$P_k'(t_{n+1})$$

近似真实导数 (y’(t_{n+1}))，并要求

$$P_k'(t_{n+1}) = f(t_{n+1},y_{n+1}).$$

由于右边包含未知的 (y_{n+1})，BDF天然是隐式方法。

对于BDF方法， 我们采用Newton后向差值多项式来进行预测。给定节点 (x_0,x_1,\ldots,x_k)，Newton插值多项式为

$$\begin{aligned} P_k(x) ={}&y[x_0] +(x-x_0)y[x_0,x_1]\\ &+(x-x_0)(x-x_1)y[x_0,x_1,x_2]\\ &+\cdots\\ &+\prod_{m=0}^{k-1}(x-x_m)y[x_0,\ldots,x_k]. \end{aligned}$$

其中 (y[x_0,\ldots,x_j]) 是 (j) 阶差商。

对于等距的后向节点，差商可以用后向差分表示。定义后向差分

$$\nabla y_{n+1}=y_{n+1}-y_n,$$ $$\nabla^2 y_{n+1} = y_{n+1}-2y_n+y_{n-1},$$

一般地，

$$\nabla^j y_{n+1} = \sum_{\ell=0}^{j} (-1)^\ell \binom{j}{\ell} y_{n+1-\ell}.$$

令

$$t=t_{n+1}+\theta h.$$

在 (t_{n+1}) 附近，Newton后向插值多项式为

$$P_k(t_{n+1}+\theta h) = \sum_{j=0}^{k} \binom{\theta+j-1}{j} \nabla^j y_{n+1}.$$

在 (\theta=0) 处对时间求导，可以得到

$$P_k'(t_{n+1}) = \frac1h \sum_{j=1}^{k}\frac1j\nabla^j y_{n+1}.$$

因此，任意 (k) 阶BDF公式为

$$\boxed{ \sum_{j=1}^{k}\frac1j\nabla^j y_{n+1} = h f(t_{n+1},y_{n+1}) }$$

这就是定步长BDF最紧凑的一般形式。

前几个BDF公式:

BDF1：后向欧拉

$$y_{n+1}-y_n=hf_{n+1}$$

BDF2

$$\frac32y_{n+1} -2y_n+\frac12y_{n-1} = hf_{n+1}$$

BDF3

$$\frac{11}{6}y_{n+1} -3y_n+\frac32y_{n-1} -\frac13y_{n-2} = hf_{n+1}$$

BDF4

$$\frac{25}{12}y_{n+1} -4y_n+3y_{n-1} -\frac43y_{n-2} +\frac14y_{n-3} = hf_{n+1}$$

BDF5

$$\frac{137}{60}y_{n+1} -5y_n+5y_{n-1} -\frac{10}{3}y_{n-2} +\frac54y_{n-3} -\frac15y_{n-4} = hf_{n+1}$$

BDF6

$$\frac{49}{20}y_{n+1} -6y_n+\frac{15}{2}y_{n-1} -\frac{20}{3}y_{n-2} +\frac{15}{4}y_{n-3} -\frac65y_{n-4} +\frac16y_{n-5} = hf_{n+1}$$

3.5 Adams和BDF的稳定性

对于Adams和BDF的稳定性分析较为复杂，笔者计划另开专题讨论…

TL;DR:

• Adams 方法：精度高、每步成本低，适合非刚性问题；但高阶时稳定区域有限，遇到刚性方程会被迫使用很小步长。
• BDF 方法：阶数通常较低、每步需要隐式求解，但在负实轴附近稳定区域很大，能有效压制快速衰减模态，因此适合刚性问题。

---

参考文献与源码

1. Linda R. Petzold, Automatic Selection of Methods for Solving Stiff and Nonstiff Systems of Ordinary Differential Equations, SIAM Journal on Scientific and Statistical Computing, 4(1), 136–148, 1983. DOI: 10.1137/0904010.
2. Alan C. Hindmarsh, ODEPACK, A Systematized Collection of ODE Solvers, in Scientific Computing, North-Holland, 1983, pp. 55–64. LLNL PDF.
3. Lawrence Livermore National Laboratory, ODEPACK 项目说明与官方软件页面.
4. ODEPACK double precision source: 主驱动源码 opkdmain.f与辅助源码 opkda1.f.
5. K. Radhakrishnan and A. C. Hindmarsh, Description and Use of LSODE, the Livermore Solver for Ordinary Differential Equations, LLNL Report UCRL-ID-113855, 1993. LLNL PDF.
6. SciPy documentation, scipy.integrate.LSODA and solve_ivp.