隐式 ODE 方法与非线性步求解

Newton 迭代与线性求解的接口 ​

对于向量状态，后向 Euler 的 Jacobian 为

第 $k$ 次 Newton 迭代理路Newton 非线性方程法Newton's method · Newton-Raphson method · Newton method for nonlinear systems在当前点解一阶线性化方程来修正非线性方程的近似解，并分析其局部二次收敛与失效边界。解

这里通常应解线性系统，而不是显式形成矩阵逆。若 $f$ 来自 PDE 空间半离散理路线法：先空间离散再时间积分Method of lines · MOL · Semi-discretization in space先把 PDE 变为空间半离散 ODE，再分析时间推进、边界输入、刚性、质量矩阵和空间与时间误差的配合。，$J_f$ 常是大型稀疏矩阵；每步效率主要取决于稀疏分解、Krylov 迭代和预条件，而不只是 Newton 外层做了几轮。

完整 Newton在每轮更新 Jacobian；简化 Newton在若干轮或若干步内复用近似矩阵。复用减少组装和分解成本，但状态或步长变化较大时，旧矩阵可能不再有效。以收敛速度和失败情况决定何时更新，比固定宣称“每步复用一定更快”更合理。[1]

不精确 Newton 允许内层线性系统只近似求解。远离根时，外层误差很大，没有必要把线性残差压得极低；接近根时，再提高线性求解精度。这样让内外误差处在匹配尺度上，而不是把某一层过度求解。

不动点迭代为什么可能失败 ​

最直接的迭代是

若 $f(t_{n+1},\cdot )$ 在整个有限维状态空间上全局 Lipschitz，且 $hL<1$，则右端是完备空间上的压缩自映射，可用Banach 不动点定理理路Banach 不动点定理Banach fixed-point theorem · Contraction mapping theorem完备空间中的统一压缩给出唯一不动点；用几何尾和证明收敛，并把后验误差与残差转成停止证书。证明收敛。若只在局部非空闭集上有 Lipschitz 估计，还必须验证该集对迭代不变；仅有 $hL<1$ 不能保证迭代留在有估计的区域。这是一个充分条件，不是隐式方程可解的必要条件。

在线性衰减 $y^{\prime }=-\lambda y$、$\lambda >0$ 上，不动点误差每轮乘以 $-h\lambda$。当 $h\lambda >1$ 时它发散，但后向 Euler 方程 $(1+h\lambda )z=y_n$ 仍然唯一可解，而且 Newton 在线性情形下一步就能求出精确根。时间格式稳定与求根迭代收敛是两个不同问题。

补充示意图

一个能完整算出来的非线性步 ​

考虑 $y^{\prime }=-y^2$，当前值 $y_n=1$，步长 $h=1/2$。后向 Euler 要求

代数方程有两个根 $-1\pm \sqrt{3}$。沿着 $h\downarrow 0$ 与旧状态连续连接的物理解支，是

另一个负根并不是这条衰减轨道的小步延续。隐式方程在非线性情况下可能存在多个根，因此“残差为零”还应与所跟踪的解支一致。良好预测值和适当步长在这里不仅影响速度，也帮助定位正确解支。

使用 Newton 法，$G^{\prime }(z)=1+z$，从 $z^{(0)}=1$ 出发：

前三次更新为

几次迭代后已经非常接近这个后向 Euler 步的精确根。可是原 ODE 在时间 $h=1/2$ 的精确解是 $1/(1+h)=2/3$。两者仍有约 $0.0654$ 的差距，这是时间离散误差，不是 Newton 还没有解准。

这个例子把三件事清楚分开了：原微分方程的解、隐式离散方程的解，以及当前迭代值。继续迭代只会消除后两者之间的差，不会把一阶时间格式变成精确积分。

残差多小才算这一步解好了 ​

设理想离散根为 $z_{\ast }$，当前值 $z$ 足够接近它，且 Jacobian 可逆。线性化给出

所以残差 $G(z)$ 与状态误差之间还隔着逆 Jacobian。残差很小而问题局部病态时，状态误差未必同样小；反过来，残差的物理单位或分量尺度很不均匀时，直接比较普通 Euclidean 范数也可能失真。

对 $m\ge 1$ 个状态分量，常见的尺度化方式是选 ${\mathrm{atol}}_i>0$、$\mathrm{rtol}\ge 0$ 并设置

绝对容差处理接近零的量，相对容差处理随状态尺度变化的量。这是一种常用形式；具体求解器可能采用当前解或预测值构造权重，应以实现约定为准。Newton 修正量、非线性残差与外层局部误差估计也不是同一个数，不能随意共用未经缩放的阈值。[3]

为什么不能永远使用固定绝对残差 ​

考虑最简单的方程 $y^{\prime }=0$。后向 Euler 本来应保持常数。若每步只算到

那么 $T/h$ 步后，累积误差为约 $T\eta /h$。即使每一步的 $\eta$ 看起来很小，保持它不变而无限缩小时间步，误差仍可能增大。

在具有稳定传播和适当逆 Jacobian 控制的情形，若一个 $p$ 阶时间方法每步额外引入 $O(h^{p+1})$ 的状态扰动，累积后可与 $O(h^p)$ 时间误差保持同阶。这里说的是状态扰动，不是任意未经缩放的残差。自适应变步长求解器会使用更细致的局部误差预算，而不是对所有问题强行套用同一个 $h^{p+1}$ 数值阈值。

求根收敛和局部时间误差合格是两道不同检查；外层误差过大时，已经收敛的候选步仍可能被拒绝。

A 稳定与 L 稳定分别解决什么 ​

用测试方程 $y^{\prime }=\lambda y$，记 $\zeta =h\lambda$。后向 Euler 的稳定函数为

当 $\mathrm{Re}\zeta \le 0$ 时，$|1-\zeta |\ge 1$，所以 $|R_{\mathrm{BE}}|\le 1$；这就是 A 稳定理路ODE 方法的绝对稳定域Absolute stability region · Region of absolute stability以标量衰减测试方程的放大因子定义一步法稳定域，并区分 A-stability、L-stability、精度与其他稳定概念。。同时当 $\zeta$ 沿负实轴趋向 $-\mathrm{\infty }$ 时，$R_{\mathrm{BE}}\to 0$，快速衰减模态会被有效压小，这体现了 L 稳定所强调的强刚性理路刚性常微分方程Stiff ordinary differential equation · Stiff ODE · Stiffness刻画快速衰减模态迫使显式方法使用远小于精度需求的步长这一方法相关现象，并说明隐式路线的能力与代价。阻尼。

梯形法的稳定函数为

它也是 A 稳定，但 $\zeta \to -\mathrm{\infty }$ 时趋于 $-1$。例如 $\zeta =-100$ 时，后向 Euler 振幅因子为 $1/101$，梯形法为 $-49/51$，而精确振幅因子是 $e^{-100}$。梯形法没有爆炸，却可能保留一个反复变号、缓慢衰减的数值分量。

这并不说明梯形法普遍较差：在所需时间尺度被充分解析时，它具有二阶精度，而后向 Euler 只有一阶。比较的重点是，当前任务需要多少阻尼、是否存在未解析的刚性瞬态，以及时间误差目标是什么。[2]

隐式方法也并非全部 A 稳定。不同 Runge–Kutta 或多步公式有不同稳定域；例如 BDF 族的高阶成员并不都具有整个左半平面的稳定性。因此选择时间方法仍需查看其具体稳定函数或稳定区域，不能把“隐式”当作一个通用稳定性证明。

一步的接受，包含两道检查 ​

一个实用隐式积分步通常先用旧解、外推或较便宜公式产生预测值，再用 Newton 或其他方法求解代数方程。只有非线性迭代达到要求，才获得一个足够准确的候选时间步。

但候选步还要通过外层局部时间误差测试。非线性方程解得很准，只说明实现了这一步公式；若时间步过大，公式本身的截断误差仍会超标。自适应算法可能因此拒绝一个 Newton 已经收敛的步，并缩小 $h$ 重算。

若 Newton 不收敛，可以更新 Jacobian、改善预测值或采用阻尼；必要时也会缩小时间步。缩步同时改变隐式方程和线性系统，往往改善局部可解性与迭代条件。接受后再根据局部误差估计调整下一步，而不是根据 Newton 轮数单独决定精度是否合格。

对具有多个物理时间尺度的问题，这样的内外两层控制是隐式积分器的重要组成部分。只给出后向 Euler 的公式，却不说明如何解根、怎样判断收敛、何时拒步，仍不足以构成一个可靠的通用算法。

更一般的隐式结构 ​

对质量矩阵系统 $M\dot{y}=F(t,y)$，后向 Euler 残差为

若 $M$ 奇异，则要按微分代数方程处理约束与一致初值。若 $M$ 依赖状态，求 Jacobian 时还必须考虑它的导数，不能直接照搬常质量矩阵公式。

微分代数系统理路微分代数方程与一致初值Differential-algebraic equation · DAE · 微分代数系统将隐式微分关系、代数约束和一致初值放在同一模型中，以可逆代数偏导块证明局部可解性，并验证后向Euler的邻近分支和约束残差。把“即时残差为零”与“确有轨道通过该状态”分开。可逆代数偏导块允许用隐函数约化，再以Schur补解释后向Euler的邻近分支；正则常系数DAE理路正则常系数微分代数系统Regular linear DAE · Constant-coefficient differential-algebraic system · Weierstrass decomposition for DAEs从可逆移位与核像分解证明正则矩阵束的微分—幂零分块，给出全部经典解和一致初值，并把外力导数损失与隐式步的误差放大写成有限证书。还给出每步矩阵可逆、错误初值却产生随缩步放大尖峰的精确例子。因此本页的隐式求解成功记录仍应另外核验约束和初始相容性。

隐式 Runge–Kutta 方法把若干阶段量作为未知数。对角隐式方法可以按阶段顺序求解；全隐式方法一般形成耦合块系统。BDF 等多步法则使用历史状态构造单个新状态方程，同时需要启动过程和变步长历史管理。这些方法的时间离散结构不同，但求解内核仍围绕残差、Jacobian、线性代数和误差预算组织。

适合的选择不是“遇到 ODE 就用最复杂的隐式法”，而是让稳定域、阻尼、精度与线性系统结构匹配。对于扩散半离散这类刚性问题，隐式方法的优势往往来自能够复用稀疏结构或预条件器，并把步长交回真正需要解析的时间尺度。

隐式方法仍遵循常微分方程时间步进框架理路常微分方程时间步进框架Time-stepping method for ODE · ODE integrator · Numerical time integration将常微分方程初值问题离散为时间网格上的数值状态与一步或多步更新，并分离流映射、误差层和退出状态。；每次 Newton 更新又归结为求解一个线性方程组理路线性方程组System of linear equations可写为矩阵方程 Ax=b 的有限个一次方程系统。。