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

[TOC]

LSODA方法中所谓的“自适应”，其实包含以下三个方面：

• 步长自适应
• 阶数自适应
• 方法族自适应

在上一篇文章中，我们已经提到LSODA使用的两种方法，但是还未提及它在实际计算的时候是如何自适应地切换的。在这篇文章中，我们将分别来讲一下它是如何做到这三者的自适应的。

4. 从定步长公式到自适应多步法

在上一篇文章的第 3 节，我们从等距节点出发，推导了定步长 Adams 与 BDF 公式。不过在实际求解的过程中，定步长所带来的缺陷是十分明显的。我们知道，实际 ODE 的时间尺度通常并不均匀。解在某些区间可能十分平滑，在另一些区间却可能出现快速过渡。如果始终使用同一个步长，就会面临以下两种问题：

• 若步长按最困难区间选取，那么在平滑区间会进行大量没有必要的计算；
• 若步长按平滑区间选取，那么进入快速变化区间后，局部误差可能超过容差，甚至发生迭代或稳定性崩溃。

从原则上说，一个可接受的步长至少同时受到精度与稳定性两方面的约束（这点我们在上一篇文章中已经举过例子了）：

$$|h| \le \min\bigl(h_{\mathrm{accuracy}},h_{\mathrm{stability}}\bigr).$$

其中 由局部截断误差决定； 则取决于方法的绝对稳定区域和当前动力学尺度。由于两者都是局部的性质，所以自适应求解器需要在每一步之后重新考虑

$$h_{n+1} = r_n h_n,$$

其中 是根据误差、稳定性、迭代收敛情况和安全因子得到的步长倍率。关于这个倍率具体是如何确定的，我们将放到后面详细讲解。

4.2 变步长给多步法带来的额外困难

我们知道 Adams 和 BDF 是多步法，当前公式依赖若干历史节点的信息。令

$$h_n=t_{n+1}-t_n,$$

一旦 不再相等，历史节点在当前步长尺度下的位置也随之改变，定步长系数一般不能原样使用。

例如，对 Adams 方法令

$$\theta=\frac{t-t_n}{h_n}, \qquad c_j=\frac{t_{n-j}-t_n}{h_n},$$

并用 Lagrange 基函数

$$L_j(\theta) = \prod_{m\ne j} \frac{\theta-c_m}{c_j-c_m}$$

插值历史导数，则相应积分权重为

$$b_j=\int_0^1L_j(\theta)\,d\theta.$$

因此一般变步长 Adams 公式中的 依赖所有相关历史步长比。

BDF 也有同样现象。以 为中心，令

$$\xi=\frac{t-t_{n+1}}{h_n}, \qquad c_j=\frac{t_{n+1-j}-t_{n+1}}{h_n},$$

其中 ，更早节点的位置由历史步长比决定。若 是这些节点上的 Lagrange 基函数（你也许还记得，在上一篇文章中，我们使用Newtown差值法来对 进行预测。但是那是定步长的情形，我们可以将向后微商写成向后差分的形式。但在这里，我们考虑了变步长，这样的情况下使用Lagrange基函数会更加合适一些），则一般变步长 BDF 可以写成

$$\boxed{ \sum_{j=0}^{q}L_j'(0)\mathbf y_{n+1-j} =h_n\mathbf f(t_{n+1},\mathbf y_{n+1}). }$$

例如令 ，变步长 BDF2 为

$$\boxed{ \frac{1+2r}{1+r}\mathbf y_{n+1} -(1+r)\mathbf y_n +\frac{r^2}{1+r}\mathbf y_{n-1} =h_n\mathbf f_{n+1}. }$$

当 时，它才退化为定步长 BDF2：

$$\frac32\mathbf y_{n+1}-2\mathbf y_n+\frac12\mathbf y_{n-1} =h\mathbf f_{n+1}.$$

上面的推导描述了将定步长的多步法推广到变步长的基本思想。可以看到，由于步长的改变，公式中的系数也会发生相应的修改。但是如果每次计算一个step都需要按照上面的方法重新计算每一个历史节点的系数，可想而知这样的效率是十分低下的。所以LSODA采用下面所介绍的局部多项式和Nordsieck向量，来避免每次改变步长后，重新处理不等距历史节点和计算整套权重。

4.3 局部多项式与 Nordsieck 向量

假设当前内部时间为 ，当前步长为 ，方法阶数为 。定义缩放导数（这里我们之所以使用约等号，是因为实际上我们并不知道严格的导数值）

$$\mathbf z_j(t_n) \approx \frac{h^j}{j!}\mathbf y^{(j)}(t_n), \qquad j=0,1,\ldots,q.$$

把它们组成 Nordsieck 向量

$$\mathbf z_n= \begin{bmatrix} \mathbf z_0\\ \mathbf z_1\\ \vdots\\ \mathbf z_q \end{bmatrix} = \begin{bmatrix} \mathbf y_n\\ h\mathbf y'_n\\ \vdots\\ \dfrac{h^q}{q!}\mathbf y_n^{(q)} \end{bmatrix}.$$

等价地，引入无量纲变量

$$\theta=\frac{t-t_n}{h},$$

局部解可以表示为

$$\mathbf Q_n(\theta) = \sum_{j=0}^{q}\mathbf z_j\theta^j.$$

对局部多项式 (Q(x))，有

$$Q(x)\approx y(t_{n+1}+hx),$$

而对 (x) 求导时

$$Q'(x)\approx h\,y'(t)=h f(t,y(t)).$$

因此：

• (Q(-1)) 表示 (y_n)；
• (Q’(-1)) 表示 (hf_n)；
• (Q’(-2)) 表示 (hf_{n-1})；
• (Q’(0)) 表示新时刻的 (hf_{n+1})。

在这个局部多项式中，我们巧妙地利用了局部多项式来储存历史信息，用修正的方式来消除每一步步长不相等的影响。

4.4 Nordsieck向量在变步长下的变换

若此时新的步长为

$$h_{\mathrm{new}}=r h,$$

定义新尺度变量

$$\eta=\frac{t-t_n}{h_{\mathrm{new}}}.$$

因为

$$\theta=r\eta,$$

所以同一条多项式可以写成

$$Q_n = \sum_{j=0}^q z_j\theta^j = \sum_{j=0}^q r^jz_j\eta^j.$$

因此

$$\widetilde z_j=r^jz_j.$$ $$(t_n,h,\mathbf z_n) \longrightarrow (t_n,h_{\mathrm{new}},D(r)\mathbf z_n)$$

其中

$$D(r)=\operatorname{diag}(1,r,r^2,\ldots,r^q).$$

注意到此时的时间中心还是 (t_n)，接下来我们就需要进行 的变换。

令

$$x=\frac{t-t_{n+1}}h,$$

则 ，所以

$$\mathbf Q_n(x+1) = \sum_{j=0}^{q}\mathbf z_j(x+1)^j.$$

把它重新按 的幂次展开，就得到新时刻的预测 Nordsieck 向量

$$\boxed{ \mathbf z_{n+1}^{(p)}=P\mathbf z_n, }$$

其中 是由二项式系数组成的上三角 Pascal 矩阵。三阶时

$$P= \begin{bmatrix} 1&1&1&1\\ 0&1&2&3\\ 0&0&1&3\\ 0&0&0&1 \end{bmatrix}.$$

至此我们就得到了在变步长情况下，多步法中历史节点的Nordsieck向量储存方式和更新方法。

4.6 Pascal平移的公共预测及其校正

预测多项式在新时刻给出

$$\mathbf y_{n+1}^{(p)}=\mathbf z_0^{(p)}, \qquad h\mathbf y_{n+1}^{\prime(p)}=\mathbf z_1^{(p)}.$$

但预测值一般不严格满足微分方程（接下来我们会看到，这里的 Pascal 平移给出的是公共预测。而之后的修正，将体现出Adams和BDF两种算法的差别）。ODE要求

$$h\mathbf y'_{n+1} = h\mathbf f(t_{n+1},\mathbf y_{n+1}).$$

因此定义新时刻的缩放导数误差

$$\boxed{ \mathbf a = h\mathbf f(t_{n+1},\mathbf y_{n+1}) -\mathbf z_1^{(p)}. }$$

若只修改一阶导数，新的多项式就会与原有的历史条件不相容（在下面一小节可以更具体地看到）。正确做法是选择一个标量校正多项式

$$l(x)=l_0+l_1x+\cdots+l_qx^q,$$

并令

$$\boxed{ \mathbf Q^{(c)}(x) = \mathbf Q^{(p)}(x)+l(x)\mathbf a. }$$

于是各阶系数统一更新为

$$\boxed{ \mathbf z_j^{(c)} = \mathbf z_j^{(p)}+l_j\mathbf a, \qquad j=0,\ldots,q. }$$

特别地，校正后的解是

$$\mathbf y_{n+1} = \mathbf z_0^{(p)}+l_0\mathbf a.$$

由于 被定义为一阶缩放导数的修正，必须要求

$$l'(0)=1,$$

也就是 。Adams 与 BDF 的差别，现在转化成了：应当选取什么样的 ，才能在修正新时刻的同时保留各自依赖的历史信息？既然我们已经把历史节点信息用局部多项式 的形式保存下来了（并且还统一了步长），那么Adams和BDF的修正也就很好推导得到了。

4.7 Adams 的校正多项式

(q) 阶隐式 Adams–Moulton 方法，本质上用以下导数值构造一个 (q-1) 次导数插值多项式：

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

在 (x) 坐标中，这些节点是

$$0,-1,-2,\ldots,-(q-1).$$

然后将这个导数插值多项式从 (t_n) 积分到 (t_{n+1})，因此还需要一个积分起点：

$$y(t_n)=y_n.$$

所以一个 (q) 阶 Adams 局部解多项式由以下 (q+1) 个条件确定：

$$\boxed{ Q(-1)=y_n }$$

以及

$$\boxed{ Q'(0)=hf_{n+1},\quad Q'(-1)=hf_n,\quad\ldots,\quad Q'(-(q-1))=hf_{n-q+2}. }$$

次数为 (q) 的多项式有 (q+1) 个系数，恰好由这些条件唯一确定。

经过 Nordsieck–Pascal 预测后，我们得到

$$Q^{(p)}(x)=Q^{(\mathrm{old})}(x+1).$$

它在新时刻给出预测值和预测导数：

$$y_{n+1}^{(p)}=Q^{(p)}(0), \qquad h{y'}_{n+1}^{(p)}={Q^{(p)}}'(0).$$

在固定步长原型中，这条预测多项式仍然携带旧的 Adams 历史条件：

$$Q^{(p)}(-1)=y_n,$$ $${Q^{(p)}}'(-j)=hf_{n+1-j}, \qquad j=1,\ldots,q-1.$$

问题发生在新时刻的导数项：

$${Q^{(p)}}'(0) \ne hf(t_{n+1},Q^{(p)}(0)).$$

我们并不能保证这两者是相等的，因为我们并没有 的信息，这是靠外推得到的。我们希望校正后的多项式满足

$${Q^{(c)}}'(0) = hf(t_{n+1},Q^{(c)}(0)).$$

于是可以定义所需的新导数改变量

$$a= hf(t_{n+1},y_{n+1}) - {Q^{(p)}}'(0).$$

但是如果我们只把一阶系数改成

$$z_1^{(c)}=z_1^{(p)}+a,$$

相当于令

$$Q^{(c)}(x)=Q^{(p)}(x)+xa.$$

但在历史节点 (x=-j) 处，

$$Q^{(c)}(-j)=Q^{(p)}(-j)-ja,$$

会破坏已有的历史插值条件。

因此，校正不能只改变一个系数，而必须沿某个完整多项式方向进行：

$$Q^{(c)}(x)=Q^{(p)}(x)+l(x)a.$$

这里的 (l(x)) 要完成两件事：

1. 在新时刻把导数增加 (a)；
2. 在历史节点不改变 Adams 所依赖的信息。

由于我们的积分起点 需要保留，要求

$$Q^{(c)}(-1)=Q^{(p)}(-1)=y_n.$$

所以必须有

$$\boxed{l(-1)=0.}$$

Adams 使用历史导数需要保持不变

$$hf_n,hf_{n-1},\ldots,hf_{n-q+2}.$$

要求校正后不改变这些导数：

$${Q^{(c)}}'(-j) = {Q^{(p)}}'(-j).$$

所以必须要求

$$\boxed{ l'(-1)=l'(-2)=\cdots=l'(-(q-1))=0. }$$

在新时刻把导数增加 (a)，要求

$${Q^{(c)}}'(0) = {Q^{(p)}}'(0)+a.$$

所以必须有

$$\boxed{l'(0)=1.}$$

利用上面对于 的约束，我们可以对其进行求解。(l(x)) 是 (q) 次多项式，因此 (l’(x)) 是 (q-1) 次多项式。我们已经知道 (l’(x)) 有 (q-1) 个零点：

$$-1,-2,\ldots,-(q-1).$$

所以它一定具有形式

$$l'(x) = C(x+1)(x+2)\cdots(x+q-1).$$

再使用归一化条件 (l’(0)=1)：

$$1 = C(1)(2)\cdots(q-1) = C(q-1)!,$$

因此

$$\boxed{ l'(x) = \frac{(x+1)(x+2)\cdots(x+q-1)} {(q-1)!}. }$$

最后，导数只能确定 (l(x)) 到一个积分常数。利用

$$l(-1)=0$$

得到

$$\boxed{ l(x) = \int_{-1}^{x} \frac{(s+1)(s+2)\cdots(s+q-1)} {(q-1)!}\,ds. }$$

这就是 Adams 校正多项式。

4.8 BDF 的校正多项式

在固定步长原型中， 阶 BDF 要保留前 个历史解值，因此要求

$$l(-1)=l(-2)=\cdots=l(-q)=0.$$

满足这些零点条件的多项式必为常数倍的

$$(x+1)(x+2)\cdots(x+q).$$

再用 归一化。由于

$$\left. \frac{d}{dx}\prod_{j=1}^{q}(x+j) \right|_{x=0} =q!H_q, \qquad H_q=\sum_{j=1}^{q}\frac1j,$$

得到

$$\boxed{ l(x) = \frac{(x+1)(x+2)\cdots(x+q)}{q!H_q}. }$$

这就是 BDF 校正多项式。

4.9 算法的等价性

看到这里或许你会有一个疑问

在这一套预测-校正算法中，我没有看到我们使用Lagrange多项式对f或者y进行预测，全程都是Q在预测，而Q的建立是根据Taylor展开得到的。那我们是如何保证我们这一套预测-校正算法最终给出的每个历史节点的权重，与我们常规定步长所使用的权重是一致的呢？

其根本原因在于，我们都在尝试用一个 阶的多项式来拟合 或者 。同一个 (q) 次多项式可以写成不同形式。

Lagrange 形式：

$$Q(x)=\sum_{j=0}^q L_j(x)y_j.$$

Newton 形式：

$$Q(x)=y[c_0]+(x-c_0)y[c_0,c_1]+\cdots.$$

Nordsieck 形式：

$$Q(x)=z_0+z_1x+\cdots+z_qx^q.$$

这三者只是同一多项式的不同坐标表示：

$$\boxed{ \text{Lagrange 节点值} \longleftrightarrow \text{Newton 差商} \longleftrightarrow \text{Nordsieck 系数}. }$$

---

5. 隐式校正方程的求解

上一节把 Adams 与 BDF 都归结为同一个非线性方程

$$\mathbf G(\mathbf a) = \mathbf a+\mathbf z_1^{(p)} -h\mathbf f\!\left( t_{n+1}, \mathbf z_0^{(p)}+l_0\mathbf a \right) =0.$$

预测值对应初始猜测 。一旦求得 ，新解和整个历史多项式便同时确定：

$$\mathbf y_{n+1}=\mathbf z_0^{(p)}+l_0\mathbf a, \qquad \mathbf z_{n+1}=\mathbf z_{n+1}^{(p)}+\mathbf l\mathbf a.$$

因此接下来我们的问题就变成了怎样求解 。

5.1 函数迭代：便宜但有收敛半径限制

直接把校正方程写成不动点形式：

$$\boxed{ \mathbf a^{(m+1)} =h\mathbf f\!\left( t_{n+1}, \mathbf z_0^{(p)}+l_0\mathbf a^{(m)} \right) -\mathbf z_1^{(p)}. }$$

设精确校正为 ，迭代误差为

$$\boldsymbol\varepsilon^{(m)} =\mathbf a^{(m)}-\mathbf a^*.$$

在解附近线性化，得到

$$\boldsymbol\varepsilon^{(m+1)} \approx h l_0J\,\boldsymbol\varepsilon^{(m)},$$

其中

$$J= \frac{\partial\mathbf f}{\partial\mathbf y} (t_{n+1},\mathbf y_{n+1}).$$

所以局部收敛要求近似满足

$$\boxed{\rho_{\mathrm{sp}}(h l_0J)<1}$$

其中 表示矩阵谱半径。

对非刚性问题，精度允许的步长通常使 足够小，局部收敛可以得到很好的满足。由于不动点形式的函数迭代只需要重复计算 ，因而十分便宜。DLSODA 的 Adams 分支采用这种策略。

5.2 Newton 与 chord 迭代：消除简单不动点限制

对 作 Newton 线性化：

$$\frac{\partial\mathbf G}{\partial\mathbf a} =I-hl_0J.$$

第 次迭代求解

$$\boxed{ \left(I-hl_0J^{(m)}\right) \Delta\mathbf a^{(m)} =-\mathbf G(\mathbf a^{(m)}), }$$

然后更新

$$\mathbf a^{(m+1)} =\mathbf a^{(m)}+\Delta\mathbf a^{(m)}.$$

严格 Newton 法在每次迭代都重新计算 并分解矩阵，代价过高。DLSODA 的 BDF 分支采用 modified Newton 或 chord 迭代：冻结一个近期的近似 Jacobian ，构造

$$P=I-hl_0\bar J,$$

对 做一次 LU 分解，然后在多次校正甚至多个时间步中复用。每次迭代只需形成残差并解

$$P\Delta\mathbf a =h\mathbf f(t_{n+1},\mathbf y) -\left(\mathbf z_1^{(p)}+\mathbf a\right).$$

这样做牺牲了严格 Newton 法的二次收敛，但把最昂贵的 Jacobian 构造和矩阵分解摊销到了多个步骤上。对于刚性问题，这通常远比把 缩小到函数迭代可以收敛更有效。

5.3 校正收敛的判据

对于上面所讲的函数迭代法，我们需要确定迭代停止的判据。换句话说，就是我们何时认为这个迭代收敛了。为此，我们先定义两个容差。对解的第 个分量，取两个非负参数：

$$\mathrm{RTOL}_i\ge0, \qquad \mathrm{ATOL}_i\ge0.$$

其中 称为相对容差； 称为绝对容差。

二者共同定义第 个分量的容差尺度

$$e_i = \mathrm{RTOL}_i|y_i|+\mathrm{ATOL}_i, \qquad e_i>0.$$

因此，对一个误差或校正向量 ，条件

$$|v_i|\le e_i$$

表示第 个分量没有超过所规定的容差尺度。当 较大时，相对项 通常占主导；当 接近零时，绝对项 防止允许误差尺度同时趋于零。

令

$$W = \operatorname{diag}\!\left( \frac1{e_1},\ldots,\frac1{e_N} \right),$$

并定义加权最大范数

$$\boxed{ \|\mathbf v\|_W = \|W\mathbf v\|_\infty = \max_{1\le i\le N} \frac{|v_i|}{e_i}. }$$

因此 是无量纲量；它表示各分量相对于自身误差预算的最大比例。采用最大范数可以保证任何一个尚未收敛的分量都不会被其他分量平均掉。

现在设第 次校正产生的增量为 ，并记其加权范数为

$$\delta_m=\|\Delta\mathbf a^{(m)}\|_W.$$

若迭代已经进入线性收敛区，可以用

$$r_m\approx\frac{\delta_m}{\delta_{m-1}}$$

估计收敛率。若 ，尚未完成的尾部校正量可粗略估为（实际上就是近似等比数列）

$$\delta_m\frac{r_m}{1-r_m}.$$

更严格地说，若从当前迭代开始存在某个 ，使后续增量满足

$$\delta_{m+k}\le r^k\delta_m, \qquad k=1,2,\ldots,$$

那么当前迭代值与真正不动点 之间的距离可以估计为

$$\begin{aligned} \|\mathbf a^*-\mathbf a^{(m)}\|_W &\le \sum_{k=1}^{\infty} \|\Delta\mathbf a^{(m+k)}\|_W\\ &\le \delta_m\sum_{k=1}^{\infty}r^k\\ &= \boxed{\delta_m\frac{r}{1-r}}. \end{aligned}$$

这相当于给出了尾部校正量增量的上限，我们需要使用这个估计量来判断迭代的进度。

设 是分配给非线性方程求解的无量纲误差预算。它通常应当取为局部离散误差预算的一小部分，使求解隐式方程产生的误差不会污染多步公式本身的截断误差。于是自然的停止条件是

$$\boxed{ \delta_m\frac{r_m}{1-r_m} \le \eta_{\mathrm{NL}}, \qquad 0\le r_m<1. }$$

若 ，迭代没有进入收敛状态；若 但尾部估计仍然过大，则需要继续迭代。

如果经过若干次迭代后仍不能满足这个条件，数学上通常有两类原因：对于简单函数迭代，可能是 太大，使不动点映射不具收缩性；对于 chord 迭代，可能是所使用的近似 Jacobian 已不能准确描述当前方程。前者需要减小步长或改用更强的隐式迭代，后者可以先更新 Jacobian；若更新后仍不收敛，也应减小步长并重新求解本步。

5.4 至此得到一个完整的单步算法

暂时不考虑误差控制和方法切换，一个时间步已经可以完整描述为：

1. 若步长改变，按 缩放 Nordsieck 各列；
2. 用 Pascal 平移得到 ；
3. 以 为初值求解校正方程；
4. 用 更新局部多项式；
5. 将 作为下一步的历史。

剩下的问题就是：校正量多大才可以接受本步？如果过大，应当把步长缩小多少？这就需要我们进行局部截断误差估计。

---

6. 局部误差估计

所谓的局部截断误差指的是，假设进入当前时间步的所有历史值都精确时， 数值方法在这一个时间步中产生的误差。这是我们判断是否要更改模型步长及其阶数的重要依据，所以我们需要想方设法来估计它。

6.1 利用校正量估计局部截断误差

对于 阶算法，预测过程对任意次数不超过 的多项式都是精确的，因此导数校正量可以写为

$$a = A_q\frac{h^{q+1}}{(q+1)!} y^{(q+1)}(t_{n+1}) +O(h^{q+2}).$$

类似地，局部截断误差也可以写为

$$\tau_{n+1} = B_q\frac{h^{q+1}}{(q+1)!} y^{(q+1)}(t_{n+1}) +O(h^{q+2}).$$

只要 ，我们就有

$$\tau_{n+1} = \frac{B_q}{A_q}a +O(h^{q+2}) \equiv \frac{a}{c_q}+O(h^{q+2}), \qquad c_q=\frac{A_q}{B_q}.$$

常数 只取决于方法族和阶数，因为 与 都只由相应的预测、校正公式决定，与当前具体解无关。

以二阶 Adams—Moulton 梯形公式为例。在固定步长且历史数据精确时，预测导数是对 的线性外推：

$$f_{n+1}^{(p)}=2f_n-f_{n-1}.$$

所以导数校正量为

$$a =h\left(f_{n+1}-2f_n+f_{n-1}\right) =h^3y^{(3)}(t_n)+O(h^4).$$

而梯形公式的局部截断误差为

$$\begin{aligned} \tau_{n+1} &= y(t_{n+1})-y(t_n) -\frac h2\left(y'(t_n)+y'(t_{n+1})\right)\\ &=-\frac{h^3}{12}y^{(3)}(t_n)+O(h^4). \end{aligned}$$

因此

$$\tau_{n+1} =-\frac1{12}a+O(h^4).$$

这就说明，我们可以通过预测导数校正量来估计局部截断误差，因此记

$$\boxed{ \widehat{\boldsymbol\tau}^{[q]}_{n+1}(h) = \frac{\mathbf a}{c_q} }$$

作为本步的局部误差估计。

第 5.3 节定义的加权最大范数不仅用于判断非线性校正是否收敛，也用于判断本步的局部截断误差是否可以接受。定义当前 阶方法的无量纲局部误差指标

$$\boxed{ \mathcal E_q(h) := \|\widehat{\boldsymbol\tau}^{[q]}_{n+1}(h)\|_W \approx \left\|\frac{\mathbf a}{c_q}\right\|_W. }$$

注意到在前面我们定义了一个 ，那是因为对于修正量计算的误差属于二级的误差，我们必须对此进行更加严格的控制。

6.3 接受或拒绝一个时间步

得到 后，就可以判断当前步长是否可接受：

$$\begin{cases} \mathcal E_q(h)\le1,&\text{接受该步},\\ \mathcal E_q(h)>1,&\text{拒绝该步，撤销预测并减小步长}. \end{cases}$$

一个时间步需要同时通过两项测试：

1. 隐式校正迭代收敛；
2. 局部误差指标不超过 1。

前者保证离散方程已经被足够准确地求解，后者保证这个离散方程本身的截断误差足够小。二者缺一不可。

---

7. 由误差模型选择步长和阶数

第 4 节解释了为什么需要改变步长，以及改变步长后怎样重标度 Nordsieck 历史量。在上一节中我们也介绍了如何通过导数校正量来估计局部误差。那么接下来就来讲讲我们是如何通过误差来选择步长和阶数的。

7.1 当前阶数下的理想步长

为了同时讨论不同阶数，记 为假设采用 阶方法、步长为 时的局部误差估计，并定义相应的无量纲误差指标

$$\boxed{ \mathcal E_p(h) := \|\widehat{\boldsymbol\tau}^{[p]}_{n+1}(h)\|_W. }$$

一个 阶方法的局部截断误差主项与 成正比。因此若把步长改为 ，并假设相邻两步内高阶导数变化不大，则

$$\mathcal E_p(rh) \approx r^{p+1}\mathcal E_p(h).$$

希望下一步满足误差边界 ，可得到理想倍率

$$r_p^{\mathrm{ideal}} = \mathcal E_p(h)^{-1/(p+1)}.$$

实际选择时通常乘一个安全因子 ，使预测误差略低于容差边界：

$$\boxed{ r_p = \sigma_p\mathcal E_p(h)^{-1/(p+1)}. }$$

对当前的 阶方法，上一节已经给出

$$\widehat{\boldsymbol\tau}^{[q]}_{n+1} \approx \frac{\mathbf a}{c_q},$$

所以

$$\boxed{ r_q = \sigma_q\mathcal E_q(h)^{-1/(q+1)}. }$$

7.2 阶次的选择

需要说明的是，阶次的改变并不是任意的。我们并不希望阶次发生太大的变化，所以我们只在 三个阶次当中进行选择。我们选择阶次的核心原则就是哪一个阶次对应的理想步长最长。为此，我们先来考察一下各阶次的理想步长。

降到 阶

当前局部多项式为

$$\mathbf Q_q(\theta) = \sum_{j=0}^{q}\mathbf z_j\theta^j.$$

若降到 阶，就必须舍弃最高阶项 。由于

$$\mathbf z_q \approx \frac{h^q}{q!}\mathbf y^{(q)}(t_n),$$

它具有 阶方法局部误差主项所需的 尺度。因此存在只依赖方法族和阶数的常数 ，使降阶方案的局部误差可估计为

$$\boxed{ \widehat{\boldsymbol\tau}^{[q-1]}_{n+1} \approx \gamma_q^{-}\mathbf z_q. }$$

相应的无量纲误差指标与候选步长倍率为

$$\mathcal E_{q-1}(h) = \|\widehat{\boldsymbol\tau}^{[q-1]}_{n+1}\|_W,$$ $$\boxed{ r_{q-1} = \sigma_{q-1}\mathcal E_{q-1}(h)^{-1/q}. }$$

指数为 ，是因为 阶方法的局部误差按 缩放。

保持 阶

保持当前阶数时，直接使用由校正量得到的估计：

$$\boxed{ \widehat{\boldsymbol\tau}^{[q]}_{n+1} \approx \frac{\mathbf a}{c_q}, \qquad r_q = \sigma_q\mathcal E_q(h)^{-1/(q+1)}. }$$

升到 阶

阶方法的局部误差主项含有

$$\frac{h^{q+2}}{(q+2)!} \mathbf y^{(q+2)}(t_{n+1}).$$

这个量没有直接包含在当前只到 阶的局部多项式中，因此必须利用最近若干步最高阶系数或校正量的变化来估计。记由这些历史变化得到的下一个高阶 Nordsieck 型估计为

$$\widehat{\mathbf z}_{q+2} \approx \frac{h^{q+2}}{(q+2)!} \mathbf y^{(q+2)}(t_{n+1}).$$

于是存在只依赖方法族和阶数的常数 ，使

$$\boxed{ \widehat{\boldsymbol\tau}^{[q+1]}_{n+1} \approx \gamma_q^{+}\widehat{\mathbf z}_{q+2}. }$$

相应地，

$$\mathcal E_{q+1}(h) = \|\widehat{\boldsymbol\tau}^{[q+1]}_{n+1}\|_W,$$ $$\boxed{ r_{q+1} = \sigma_{q+1}\mathcal E_{q+1}(h)^{-1/(q+2)}. }$$

升阶估计只有在已经积累了足够的平滑历史信息时才可靠。如果尚不能得到可信的 ，这一轮就不应考虑升阶。

比较三个候选方案

最终比较

$$r_{q-1},\qquad r_q,\qquad r_{q+1},$$

选择预计允许最大下一步长的阶数。

一般情况下，升阶其实并不一定是最优的选择。提高阶数的好处是误差随 更快下降，但同时也会带来：

• 更长的历史建立时间；
• 对解的光滑性更高的要求；
• 更小或形状更不利的稳定区域；
• 更多寄生模态与舍入误差影响；
• 阶数和步长改变后更复杂的历史调整。

因此当解刚经历快速暂态、右端不够光滑或步长频繁变化时，低阶方法可能允许更大的可靠步长。

7.4 Adams 的稳定性步长上限约束

上面的 都只来自局部误差估计。对 Adams 分支，还必须检查函数迭代和 Adams 稳定域是否允许这样的步长。

设从校正迭代得到的局部 Jacobian/Lipschitz 尺度估计为 （集体估计方法可以参考8.1小节）。对每个 Adams 阶数，我们都定义一个稳定性界 ，并要求候选新步长满足近似条件

$$|h|r\widehat L\lesssim s_q.$$

因此实际允许的倍率是

$$r_{\mathrm{Adams}} = \min\left( r_{\mathrm{accuracy}}, \frac{s_q}{|h|\widehat L} \right).$$

如果第二项更小，就说明 Adams 的步长不是由误差容限决定，而是被稳定性或函数迭代收敛性限制。后面我们会看到，这正是“当前区域表现出刚性”的可计算信号。

你可能会好奇 BDF 为什么没有相应的约束，其实 BDF 也有相应的稳定性问题，但是相比于 Adams，BDF 的稳定性要求宽泛很多，这就涉及到了两者的稳定性分析。分析的具体内容我计划另写一篇文章，敬请期待！

7.5 迟滞与失败恢复

误差和高阶导数估计本身带有噪声。如果每一步都改变 和 ，反而会扰乱历史并增加拒步次数。因此实际算法加入迟滞和失败恢复（这里所谓的失败，即可能是收敛性没有达到要求，也有可能是局部截断误差没有达到要求）：

• 建议步长增加不足约 时，通常保持原步长；
• 改变步长后，延迟若干步再考虑变阶；
• 失败后限制下一次允许的放大倍率；
• 连续误差失败时更激进地缩短步长；
• 多次失败后降到一阶并重新计算导数历史。

---

8. 从稳定性受限到自动切换方法

LSODA自动切换方法的原则是：

如果 Adams 因稳定性只能使用远小于精度所允许的步长，并且 BDF 预计能以足够大的步长优势抵消矩阵运算成本，那么当前区域值得切换到 BDF。

为此，我们需要估计 Adams 方法的稳定性限制，以及两种方法的估计步长。

8.1 从函数迭代收敛率估计局部动力学尺度

Adams 函数迭代在线性化后满足

$$\boldsymbol\varepsilon^{(m+1)} \approx h l_0J\boldsymbol\varepsilon^{(m)}.$$

记连续校正增量的加权范数为

$$\delta_m=\|\Delta\mathbf a^{(m)}\|_W,$$

则比值

$$r_m \approx \frac{\delta_m}{\delta_{m-1}}$$

包含了 的局部放大信息。因此可以构造

$$\widehat L \approx \frac{r_m}{|h l_0|},$$

从7.4小节我们知道，Adams 步长的选取收到稳定性的限制，于是我们可以据此来判断 Adams 的候选步长是否会被稳定性压制。这点是十分重要的，因为在比较两者的步长时，Adams的估计步长需要考虑稳定性的影响（下面我们还会再特意强调这一点）。

这一估计几乎是免费的，因为校正迭代本来就需要计算 。

8.2 两种方法怎样在同一尺度上比较

定义共同的局部高阶导数尺度

$$\boldsymbol\chi_{q+1}(h) := \frac{h^{q+1}}{(q+1)!} \mathbf y^{(q+1)}(t_{n+1}).$$

所谓的“共同”，指的是不论对于 Adams 方法还是 BDF 方法，其误差表达式中都含有这样的因子。

我们可以利用当前方法的校正量给出它的估计：

$$\widehat{\boldsymbol\chi}_{q+1}(h) = \frac{\mathbf a_M}{A_{M,q}}.$$

再将其分别代入两种方法的误差主项，就得到

$$\widehat{\boldsymbol\tau}_A^{[q]}(h) = B_{A,q}\widehat{\boldsymbol\chi}_{q+1}(h), \qquad \widehat{\boldsymbol\tau}_B^{[q]}(h) = B_{B,q}\widehat{\boldsymbol\chi}_{q+1}(h).$$

于是，两种方法在同一误差权重下的无量纲局部误差为

$$\mathcal E_{A,q}(h) = \left\| \widehat{\boldsymbol\tau}_A^{[q]}(h) \right\|_W, \qquad \mathcal E_{B,q}(h) = \left\| \widehat{\boldsymbol\tau}_B^{[q]}(h) \right\|_W.$$

实际上，我们可以直接利用

$$\widehat{\boldsymbol\tau}_B^{[q]}(h) \approx \frac{B_{B,q}}{B_{A,q}} \widehat{\boldsymbol\tau}_A^{[q]}(h).$$

来估计另一方的误差。

得到误差估计后，利用 阶方法的误差缩放关系

$$\mathcal E_{M,q}(rh) \approx r^{q+1}\mathcal E_{M,q}(h),$$

可分别得到两种方法由精度允许的候选步长倍率

$$r_{A,\mathrm{acc}} = \sigma_A\mathcal E_{A,q}(h)^{-1/(q+1)}, \qquad r_{B,\mathrm{acc}} = \sigma_B\mathcal E_{B,q}(h)^{-1/(q+1)},$$

其中 是安全因子。Adams 的实际候选倍率还要受到第 7.4 节的稳定性或函数迭代收敛性限制，即

$$r_A = \min\left(r_{A,\mathrm{acc}},r_{A,\mathrm{stab}}\right), \qquad r_B\approx r_{B,\mathrm{acc}}.$$

需要注意的是，这里的 只是切换前的预测。真正切换到 BDF 后，需要重新通过非线性校正收敛和局部误差检验。

如果两种候选方法的阶数不同，便不能只使用 这一个比值。此时还需要从 Nordsieck 历史和当前校正量中，分别估计候选阶数所需的

$$\boldsymbol\chi_{p+1}(h) = \frac{h^{p+1}}{(p+1)!} \mathbf y^{(p+1)}(t_{n+1}),$$

再乘以各自的误差常数。从这里我们也可以看出 Nordsieck 向量表示的优越性。

8.3 切换判据的不对称性

Adams 每步通常只需要少量右端函数计算；BDF 还需要 Jacobian、矩阵分解和线性系统求解。所以同样的步长并不意味着同样的成本。所以并不意味着只要 BDF 的步长比 Adams 的大，我们就进行方法切换。

实际上，LSODA 默认使用如下迟滞策略：

Adams 到 BDF

只有当

$$\boxed{r_B\ge5r_A}$$

时才切换。BDF 必须表现出显著的步长优势，才值得支付刚性分支的额外成本。若当前 Adams 阶数高于 BDF 的最高阶 5，则不做常规切换比较。另外，如果误差或收敛率已经接近舍入噪声，普通倍率估计可能失真。此时只有在 Adams 步长已经明确被稳定性上限截断时，才把它作为切换证据。

BDF 到 Adams

反方向的默认条件近似为

$$\boxed{r_A\ge r_B.}$$

因为 Adams 单步更便宜，只要它能够采用不小于 BDF 的步长，就可能值得切回。同样的，切换前仍需确认误差估计没有被舍入噪声淹没，并检查 Adams 的稳定性上限。

除此之外，LSODA的方法切换存在一个“保护期”。由于方法切换会改变：

• 校正多项式 ；
• 误差常数；
• 最高允许阶数；
• 非线性迭代方式；
• 是否需要 Jacobian 和矩阵工作区。

切换后立即根据尚未稳定的新历史再次判断，容易产生 Adams/BDF 来回抖动。所以LSODA 在初始化和每次切换后设置约 20 个成功步的保护计数，先积累可靠的历史数据，再重新比较两种方法。

到这里，DLSODA 的核心算法已经从基本原理构造完成。下一篇文章我们会关注程序上的实现（希望会有下一篇文章吧），因为我们知道将算法原理转化为程序实现，并不是一个简单的映射，其中还有非常多的小巧思。