数值积分是物理仿真的心脏。在具身智能的语境下，积分器的选择直接影响着仿真的稳定性、精度和速度——这三者共同决定了强化学习训练的效率和最终策略的质量。本章将深入探讨各类积分方法的数学原理、稳定性特征及其在机器人仿真中的实际应用。我们将特别关注积分器如何影响接触力的计算、约束的满足以及长时间仿真的能量守恒性。通过本章学习，您将掌握如何根据具体任务需求选择和调优积分器，以及如何在训练速度和物理真实性之间找到最佳平衡点。
积分器的核心任务是求解常微分方程组 $\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, t)$。在刚体动力学中，状态向量 $\mathbf{x}$ 通常包含位置、姿态、线速度和角速度，而 $\mathbf{f}$ 则编码了牛顿-欧拉方程。显式方法直接使用当前状态计算导数，而隐式方法则需要求解涉及未来状态的方程。
最简单的显式积分器，更新规则为：
\[\mathbf{x}_{n+1} = \mathbf{x}_n + h\mathbf{f}(\mathbf{x}_n, t_n)\]其中 $h$ 是时间步长。前向欧拉法的局部截断误差为 $O(h^2)$，全局误差为 $O(h)$。
对于简谐振子 $\ddot{x} + \omega^2 x = 0$，将其改写为一阶系统： \(\begin{bmatrix} \dot{x} \\ \dot{v} \end{bmatrix} = \begin{bmatrix} v \\ -\omega^2 x \end{bmatrix}\)
应用前向欧拉得到传递矩阵： \(\begin{bmatrix} x_{n+1} \\ v_{n+1} \end{bmatrix} = \begin{bmatrix} 1 & h \\ -h\omega^2 & 1 \end{bmatrix} \begin{bmatrix} x_n \\ v_n \end{bmatrix}\)
| 传递矩阵的特征值为 $\lambda_{1,2} = 1 \pm ih\omega$，其模为 $ | \lambda | = \sqrt{1 + h^2\omega^2} > 1$，表明能量会随时间增长——这是前向欧拉的典型特征。 |
后向欧拉使用未来状态的导数：
\[\mathbf{x}_{n+1} = \mathbf{x}_n + h\mathbf{f}(\mathbf{x}_{n+1}, t_{n+1})\]这需要求解非线性方程组，通常使用牛顿-拉夫逊迭代： \(\mathbf{x}_{n+1}^{(k+1)} = \mathbf{x}_{n+1}^{(k)} - \left(\mathbf{I} - h\frac{\partial \mathbf{f}}{\partial \mathbf{x}}\right)^{-1} \left(\mathbf{x}_{n+1}^{(k)} - \mathbf{x}_n - h\mathbf{f}(\mathbf{x}_{n+1}^{(k)}, t_{n+1})\right)\)
对于简谐振子，后向欧拉的传递矩阵为： \(\begin{bmatrix} x_{n+1} \\ v_{n+1} \end{bmatrix} = \frac{1}{1 + h^2\omega^2} \begin{bmatrix} 1 & h \\ -h\omega^2 & 1 \end{bmatrix} \begin{bmatrix} x_n \\ v_n \end{bmatrix}\)
| 特征值模为 $ | \lambda | = \frac{1}{\sqrt{1 + h^2\omega^2}} < 1$，表明能量会衰减——后向欧拉具有数值阻尼效应。 |
半隐式欧拉（也称为辛欧拉或半向后欧拉）结合了显式和隐式的特点：
\(\mathbf{v}_{n+1} = \mathbf{v}_n + h\mathbf{a}(\mathbf{x}_n, \mathbf{v}_n, t_n)\) \(\mathbf{x}_{n+1} = \mathbf{x}_n + h\mathbf{v}_{n+1}\)
这种方法先显式更新速度，再用新速度更新位置。对于分离型哈密顿系统 $H = T(\mathbf{p}) + V(\mathbf{q})$，半隐式欧拉是辛积分器的最简单例子。
考虑测试方程 $\dot{y} = \lambda y$，其中 $\lambda \in \mathbb{C}$。积分器的稳定性区域定义为复平面上使得数值解有界的 $h\lambda$ 值集合。
| 前向欧拉的稳定性条件：$ | 1 + h\lambda | \leq 1$，对应以 $(-1, 0)$ 为圆心、半径为 1 的圆盘。 |
| 后向欧拉的稳定性条件：$ | 1 - h\lambda | ^{-1} \leq 1$，对应整个左半平面的补集——这就是 A-稳定性。 |
对于刚性系统（特征值相差悬殊），隐式方法的优势明显。例如，弹簧-阻尼系统： \(\ddot{x} + 2\zeta\omega_n\dot{x} + \omega_n^2 x = 0\)
当 $\zeta \gg 1$ 时，系统刚性。显式方法需要极小的步长以保持稳定，而隐式方法可以使用较大步长。
哈密顿系统的相流保持相空间体积（刘维尔定理）和辛结构。辛积分器是保持这种几何结构的数值方法，对长时间仿真至关重要。
哈密顿系统的演化方程： \(\dot{\mathbf{q}} = \frac{\partial H}{\partial \mathbf{p}}, \quad \dot{\mathbf{p}} = -\frac{\partial H}{\partial \mathbf{q}}\)
定义辛矩阵 $\mathbf{J} = \begin{bmatrix} \mathbf{0} & \mathbf{I} \ -\mathbf{I} & \mathbf{0} \end{bmatrix}$，则哈密顿方程可写为： \(\dot{\mathbf{z}} = \mathbf{J}\nabla H(\mathbf{z})\)
其中 $\mathbf{z} = [\mathbf{q}^T, \mathbf{p}^T]^T$。
相流 $\boldsymbol{\phi}_t$ 的雅可比矩阵 $\mathbf{M} = \frac{\partial \boldsymbol{\phi}_t}{\partial \mathbf{z}}$ 满足辛条件： \(\mathbf{M}^T\mathbf{J}\mathbf{M} = \mathbf{J}\)
这保证了相空间中的”面积”（更准确地说，辛形式）守恒。
Verlet积分器是分子动力学的经典方法，有多种等价形式。
位置Verlet： \(\mathbf{x}_{n+1} = 2\mathbf{x}_n - \mathbf{x}_{n-1} + h^2\mathbf{a}_n\)
这个公式不显式包含速度，速度通过差分近似： \(\mathbf{v}_n = \frac{\mathbf{x}_{n+1} - \mathbf{x}_{n-1}}{2h}\)
速度Verlet： \(\mathbf{x}_{n+1} = \mathbf{x}_n + h\mathbf{v}_n + \frac{h^2}{2}\mathbf{a}_n\) \(\mathbf{v}_{n+1} = \mathbf{v}_n + \frac{h}{2}(\mathbf{a}_n + \mathbf{a}_{n+1})\)
速度Verlet是自启动的，且速度和位置具有相同的精度 $O(h^2)$。
Leapfrog方法在半整数时间点计算速度： \(\mathbf{v}_{n+1/2} = \mathbf{v}_{n-1/2} + h\mathbf{a}_n\) \(\mathbf{x}_{n+1} = \mathbf{x}_n + h\mathbf{v}_{n+1/2}\)
这种交错格式自然保持辛结构。对于分离型哈密顿量 $H = T(\mathbf{p}) + V(\mathbf{q})$，Leapfrog等价于以下算子分裂： \(e^{h\mathcal{L}} \approx e^{h\mathcal{L}_V/2} e^{h\mathcal{L}_T} e^{h\mathcal{L}_V/2}\)
其中 $\mathcal{L}_T$ 和 $\mathcal{L}_V$ 分别是动能和势能部分的刘维尔算子。
虽然辛积分器不精确保持能量，但能量误差有界。对于Verlet方法，存在修正哈密顿量 $\tilde{H}$： \(\tilde{H} = H + h^2H_2 + h^4H_4 + \cdots\)
使得数值解精确保持 $\tilde{H}$。这解释了为什么辛积分器的能量误差不会累积增长。
具体地，对于简谐振子，Verlet方法的能量误差为： \(\Delta E = E_{n+1} - E_n = O(h^4)\)
而长时间平均能量误差仅为 $O(h^2)$。
相空间轨迹比较
精确解 前向欧拉 Verlet
p p p
↑ ↑ ↑
│ ╱─╲ │ ╱───╲ │ ╱─╲
│ │ │ │ ╱ ╲ │ │ │
────┼───┴───┴──→ ────┼──┴───────┴─→────┼───┴───┴──→
│ q │ q │ q
│ 封闭轨道 │ 螺旋发散 │ 近似封闭
固定步长方法在处理多尺度问题时效率低下。变步长方法根据局部误差估计动态调整步长，在保证精度的同时提高效率。
Richardson外推利用不同步长的结果估计误差。设 $y(h)$ 是步长为 $h$ 的数值解，$y_{exact}$ 是精确解，若方法具有 $p$ 阶精度： \(y(h) = y_{exact} + Ch^p + O(h^{p+1})\)
使用步长 $h$ 和 $h/2$ 计算两次，误差估计为： \(\epsilon \approx \frac{|y(h/2) - y(h)|}{2^p - 1}\)
经典四阶Runge-Kutta (RK4)： \(\mathbf{k}_1 = h\mathbf{f}(t_n, \mathbf{x}_n)\) \(\mathbf{k}_2 = h\mathbf{f}(t_n + h/2, \mathbf{x}_n + \mathbf{k}_1/2)\) \(\mathbf{k}_3 = h\mathbf{f}(t_n + h/2, \mathbf{x}_n + \mathbf{k}_2/2)\) \(\mathbf{k}_4 = h\mathbf{f}(t_n + h, \mathbf{x}_n + \mathbf{k}_3)\) \(\mathbf{x}_{n+1} = \mathbf{x}_n + \frac{1}{6}(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4)\)
嵌入式Runge-Kutta方法（如Dormand-Prince）同时计算不同阶数的解，用于误差估计： \(\mathbf{x}_{n+1}^{(5)} = \mathbf{x}_n + \sum_{i=1}^7 b_i^{(5)}\mathbf{k}_i\) \(\mathbf{x}_{n+1}^{(4)} = \mathbf{x}_n + \sum_{i=1}^7 b_i^{(4)}\mathbf{k}_i\) \(\epsilon = ||\mathbf{x}_{n+1}^{(5)} - \mathbf{x}_{n+1}^{(4)}||\)
步长调整策略基于误差估计 $\epsilon$ 和容差 $\tau$： \(h_{new} = h_{old} \cdot \min\left(\alpha_{max}, \max\left(\alpha_{min}, \beta\left(\frac{\tau}{\epsilon}\right)^{1/(p+1)}\right)\right)\)
其中 $\alpha_{min}$、$\alpha_{max}$ 限制步长变化率，$\beta \approx 0.8-0.9$ 是安全系数。
对于刚体系统，还需考虑：
对于具有不同时间尺度的系统组件，可使用多速率方法：
\(\dot{\mathbf{x}}_{fast} = \mathbf{f}_{fast}(\mathbf{x}_{fast}, \mathbf{x}_{slow})\) \(\dot{\mathbf{x}}_{slow} = \mathbf{f}_{slow}(\mathbf{x}_{fast}, \mathbf{x}_{slow})\)
快变量用小步长 $h_{fast}$，慢变量用大步长 $h_{slow} = n \cdot h_{fast}$。
数值稳定性决定了仿真是否会”爆炸”。对于物理仿真，稳定性约束主要来自三个方面：显式积分的CFL条件、刚性系统的特征值分布、以及接触力的高频振荡。
考虑线性化系统 $\dot{\mathbf{x}} = \mathbf{A}\mathbf{x}$，其解为 $\mathbf{x}(t) = e^{\mathbf{A}t}\mathbf{x}_0$。设 $\mathbf{A}$ 的特征值为 $\lambda_i$，则稳定性要求 $\text{Re}(\lambda_i) \leq 0$。
对于显式欧拉，稳定性条件为： \(|1 + h\lambda_i| \leq 1, \quad \forall i\)
这要求 $h\lambda_i$ 位于以 $(-1, 0)$ 为圆心的单位圆内。对于实特征值，条件简化为： \(h \leq \frac{2}{|\lambda_{max}|}\)
对于质量-弹簧系统 $m\ddot{x} + kx = 0$，最大特征值 $\lambda_{max} = \omega = \sqrt{k/m}$，因此： \(h \leq \frac{2}{\omega} = 2\sqrt{\frac{m}{k}}\)
Courant-Friedrichs-Lewy (CFL) 条件源于波动方程的数值求解。对于波速 $c$ 和空间离散化步长 $\Delta x$，CFL数定义为： \(\text{CFL} = \frac{c \cdot h}{\Delta x}\)
稳定性要求 $\text{CFL} \leq 1$，即信息传播的数值域必须包含物理域。
在弹性体仿真中，波速由材料属性决定：
其中 $E$ 是杨氏模量，$\nu$ 是泊松比，$\rho$ 是密度，$G$ 是剪切模量。
对于显式时间积分的有限元方法： \(h \leq \alpha \frac{\Delta x_{min}}{c_{max}}\)
其中 $\alpha \approx 0.5-0.9$ 是安全系数，$\Delta x_{min}$ 是最小单元尺寸。
| 刚性系统的特征值跨度很大：$\kappa(\mathbf{A}) = | \lambda_{max} | / | \lambda_{min} | \gg 1$。常见例子： |
1. 刚性关节与柔性关节混合： 机器人系统中，某些关节刚度远大于其他关节。设关节刚度为 $k_i$，惯量为 $I_i$： \(\omega_i = \sqrt{k_i/I_i}\)
最大频率 $\omega_{max}$ 限制了全局步长。
2. 接触刚度： Hertz接触模型的接触刚度： \(k_{contact} = \frac{4}{3}E^*\sqrt{R^*\delta}\)
其中 $\delta$ 是穿透深度。接触刚度随穿透深度变化，导致变刚性问题。
3. 数值阻尼策略： 引入人工阻尼以改善条件数： \(\mathbf{M}\ddot{\mathbf{x}} + (\mathbf{C} + \alpha\mathbf{M} + \beta\mathbf{K})\dot{\mathbf{x}} + \mathbf{K}\mathbf{x} = \mathbf{f}\)
Rayleigh阻尼参数 $\alpha$、$\beta$ 的选择需要平衡稳定性和精度。
对于部分刚性系统，可将方程分解为刚性和非刚性部分： \(\dot{\mathbf{x}} = \mathbf{f}_{stiff}(\mathbf{x}) + \mathbf{f}_{non-stiff}(\mathbf{x})\)
IMEX方法对刚性部分使用隐式积分，对非刚性部分使用显式积分： \(\mathbf{x}_{n+1} = \mathbf{x}_n + h\mathbf{f}_{non-stiff}(\mathbf{x}_n) + h\mathbf{f}_{stiff}(\mathbf{x}_{n+1})\)
这在保持稳定性的同时减少了计算开销。
在强化学习训练中，实时性要求固定步长。稳定性保证策略包括：
1. 保守步长选择： \(h = \min\left(h_{target}, \gamma \cdot h_{CFL}, \frac{2}{\omega_{max}}\right)\)
其中 $\gamma \approx 0.5-0.8$ 是安全系数。
2. 子步进(Substepping)： 当检测到高频事件时，在固定的外部步长内执行多个小步长：
for external_step in range(N_external):
n_substeps = compute_required_substeps()
h_sub = h_external / n_substeps
for i in range(n_substeps):
integrate(h_sub)
3. 位置修正： 对于约束违背和穿透，使用后处理修正： \(\mathbf{x}_{corrected} = \mathbf{x}_{predicted} + \Delta\mathbf{x}_{constraint}\)
在强化学习环境中，积分器的选择直接影响策略学习的效率和质量。不当的积分器可能导致训练不稳定、样本效率低下，甚至学习到物理不真实的策略。
奖励函数通常包含位置、速度和加速度项。数值误差会导致奖励信号噪声：
能量漂移的影响： 非辛积分器的能量漂移会产生虚假的奖励信号。例如，在平衡任务中： \(r = -\alpha||\mathbf{x} - \mathbf{x}_{target}||^2 - \beta||\dot{\mathbf{x}}||^2 - \gamma||E - E_0||\)
能量项 $E$ 的数值漂移会误导策略学习。
相位误差的累积： 周期性任务（如步态）中，相位误差导致奖励时序错位： \(\phi_{numerical} = \phi_{true} + \epsilon(t)\)
长时间仿真后，$\epsilon(t)$ 可能超过 $2\pi$，完全破坏周期性。
强化学习需要大量样本，仿真速度至关重要。典型的权衡考虑：
1. 显式 vs 隐式：
2. 固定 vs 变步长：
3. 精度阶数选择：
某些积分器会引入系统性偏差，影响策略的真实性：
数值阻尼： 隐式欧拉的人工阻尼可能导致策略过度依赖这种”免费”的能量耗散： \(E_{n+1} = (1 - \xi h)E_n\)
其中 $\xi > 0$ 是数值阻尼率。
数值刚度： 某些积分器会改变系统的有效刚度： \(k_{effective} = k_{physical}(1 + \eta h^2)\)
这会影响接触力的计算和控制策略。
根据任务类型选择积分器：
操作任务（抓取、装配）：
运动任务（行走、跑步）：
软体操作：
Loup Verlet在1967年发表的论文”Computer ‘Experiments’ on Classical Fluids”标志着分子动力学仿真的重大突破。他提出的积分方法不仅计算效率高，而且具有优异的长时间稳定性。
1960年代，计算机刚刚具备模拟数百个粒子的能力。主要挑战包括：
之前的方法（如Runge-Kutta）虽然局部精度高，但长时间仿真会出现能量漂移。
Verlet的关键洞察是利用时间反演对称性。他注意到牛顿方程： \(m\ddot{\mathbf{r}} = \mathbf{F}(\mathbf{r})\)
在时间反演 $t \to -t$ 下不变。这启发他构造保持此对称性的数值格式： \(\mathbf{r}(t+h) + \mathbf{r}(t-h) = 2\mathbf{r}(t) + h^2\mathbf{a}(t) + O(h^4)\)
这个看似简单的公式隐含了深刻的几何性质——它是辛积分器。
Verlet算法的影响远超分子动力学：
其成功揭示了一个重要原则：保持物理系统的几何结构比提高局部精度更重要。
现代几何积分理论将物理系统视为流形上的动力系统，积分器应保持系统的内在几何结构。
刚体的位形空间是特殊欧几里得群 $SE(3) = SO(3) \ltimes \mathbb{R}^3$。传统方法将旋转矩阵视为 $\mathbb{R}^{9}$ 中的点，导致数值漂移： \(\mathbf{R}^T\mathbf{R} = \mathbf{I} + \epsilon(t)\)
李群积分器在群流形上演化： \(\mathbf{R}_{n+1} = \mathbf{R}_n \exp(h\boldsymbol{\omega}_n^{\times})\)
其中 $\boldsymbol{\omega}^{\times}$ 是角速度的叉乘矩阵，$\exp$ 是矩阵指数（Rodrigues公式）： \(\exp(\boldsymbol{\theta}^{\times}) = \mathbf{I} + \frac{\sin||\boldsymbol{\theta}||}{||\boldsymbol{\theta}||}\boldsymbol{\theta}^{\times} + \frac{1-\cos||\boldsymbol{\theta}||}{||\boldsymbol{\theta}||^2}(\boldsymbol{\theta}^{\times})^2\)
变分积分器从离散作用量原理出发： \(S_d = \sum_{n=0}^{N-1} L_d(\mathbf{q}_n, \mathbf{q}_{n+1}, h)\)
其中 $L_d$ 是离散拉格朗日量。离散Euler-Lagrange方程： \(D_2L_d(\mathbf{q}_{n-1}, \mathbf{q}_n, h) + D_1L_d(\mathbf{q}_n, \mathbf{q}_{n+1}, h) = 0\)
对于自由粒子 $L = \frac{1}{2}m\dot{\mathbf{q}}^2$，选择： \(L_d(\mathbf{q}_n, \mathbf{q}_{n+1}, h) = \frac{m}{2h}||\mathbf{q}_{n+1} - \mathbf{q}_n||^2\)
恢复了Verlet积分器。
对于约束系统 $\boldsymbol{\phi}(\mathbf{q}) = 0$，SHAKE和RATTLE算法保持约束流形：
SHAKE算法：
| 投影到约束流形：$\mathbf{q}{n+1} = \arg\min{\boldsymbol{\phi}(\mathbf{q})=0} | \mathbf{q} - \tilde{\mathbf{q}}_{n+1} | $ |
投影步骤通过迭代求解： \(\mathbf{q}_{n+1}^{(k+1)} = \mathbf{q}_{n+1}^{(k)} - \mathbf{G}^T(\mathbf{G}\mathbf{M}^{-1}\mathbf{G}^T)^{-1}\boldsymbol{\phi}(\mathbf{q}_{n+1}^{(k)})\)
其中 $\mathbf{G} = \partial\boldsymbol{\phi}/\partial\mathbf{q}$ 是约束雅可比。
现代观点将各种保结构性质统一在几何框架下：
这些性质通过离散微分形式和离散外微分理论严格刻画。
本章系统介绍了物理仿真中的数值积分方法，从基础的欧拉法到高级的几何积分器。核心要点包括：
基本概念与方法：
辛积分与能量守恒：
稳定性分析：
实际应用考虑：
高级技术：
关键公式汇总：
前向欧拉：$\mathbf{x}_{n+1} = \mathbf{x}_n + h\mathbf{f}(\mathbf{x}_n)$
后向欧拉：$\mathbf{x}{n+1} = \mathbf{x}_n + h\mathbf{f}(\mathbf{x}{n+1})$
Verlet：$\mathbf{x}{n+1} = 2\mathbf{x}_n - \mathbf{x}{n-1} + h^2\mathbf{a}_n$
CFL条件：$h \leq \alpha \frac{\Delta x_{min}}{c_{max}}$
李群更新：$\mathbf{R}_{n+1} = \mathbf{R}_n \exp(h\boldsymbol{\omega}_n^{\times})$
选择积分器时，需在精度、稳定性和计算效率间权衡。理解各方法的数学原理和几何意义，是设计高质量物理仿真系统的基础。
练习 2.1：证明前向欧拉方法应用于简谐振子 $\ddot{x} + \omega^2 x = 0$ 时，能量会指数增长。计算能量增长率与步长 $h$ 的关系。
提示：将二阶方程改写为一阶系统，分析传递矩阵的特征值。
练习 2.2：推导速度Verlet方法的局部截断误差。证明位置误差为 $O(h^3)$，速度误差也为 $O(h^3)$。
提示：使用泰勒展开，注意加速度在半步时刻的处理。
练习 2.3：对于刚性ODE $\dot{y} = -1000(y - \cos t) - \sin t$，比较前向欧拉和后向欧拉在不同步长下的表现。确定各方法的最大稳定步长。
提示：这是具有大负特征值的线性系统，分析稳定性区域。
练习 2.4：实现并比较RK4和Dormand-Prince自适应方法求解洛伦兹系统。记录达到给定精度所需的函数评估次数。
提示：洛伦兹系统在某些参数下呈现混沌行为，是测试积分器的好例子。
练习 2.5：设计一个保持能量和动量同时守恒的积分器用于N体问题。分析其计算复杂度和实际性能。
提示：考虑能量-动量方法或离散梯度方法。
练习 2.6：分析强化学习训练中，积分器的数值阻尼如何影响最优策略。设计实验量化这种影响。
提示：考虑摆锤摆起任务，比较不同积分器下学到的策略。
练习 2.7：推导并实现$SO(3)$上的Cayley变换积分器，比较与矩阵指数方法的精度和效率。
提示：Cayley变换 $(I-A)^{-1}(I+A)$ 将反对称矩阵映射到正交矩阵。
练习 2.8：设计一个混合积分策略，对刚体使用辛积分，对约束使用隐式方法，对接触使用显式方法。分析稳定性条件。
提示：使用算子分裂和不同时间尺度。
陷阱：使用固定的”魔法数字”步长（如always 0.01）
h = 0.01 # 错误：忽略系统特性
正确做法：基于系统特性计算步长
h = min(0.01, # 目标步长
0.5 * sqrt(m_min/k_max), # 弹簧稳定性
0.8 * dx_min/c_wave) # CFL条件
陷阱：直接积分旋转矩阵的9个元素
R = R + h * omega_cross * R # 错误：破坏正交性
正确做法：使用指数映射或重正交化
R = R @ exp(h * omega_cross) # 保持SO(3)
# 或定期重正交化
if frame % 100 == 0:
R = orthogonalize(R)
陷阱：期望辛积分器精确保持能量
assert abs(E_new - E_old) < 1e-10 # 错误：过严格
正确做法：监控长期能量有界性
# 辛积分器：能量振荡但有界
assert abs(E_avg - E_0) < tolerance * time
陷阱：固定迭代次数，不检查收敛
for i in range(5): # 错误：可能未收敛
x = implicit_step(x)
正确做法：自适应迭代with收敛检查
for i in range(max_iter):
x_new = implicit_step(x)
if norm(x_new - x) < tol:
break
x = x_new
else:
# 未收敛：减小步长或切换方法
h *= 0.5
陷阱：对所有系统使用显式方法
# 弹簧刚度 k = 10000
v += h * (-k * x / m) # 错误：需要极小步长
正确做法：识别刚性并使用合适方法
stiffness_ratio = k_max / k_min
if stiffness_ratio > 100:
use_implicit_method()
陷阱：改变步长时不更新相关状态
h = compute_new_stepsize()
x += h * v # 错误：v可能是用旧步长计算的
正确做法：同步更新所有依赖步长的量
h_new = compute_new_stepsize()
# 重新计算速度或使用插值
v = interpolate_velocity(h_old, h_new, v_old)
x += h_new * v