4. 事件检测的数值方法初步

在上一节中，我们把“事件”定义为某个标量事件函数 g(t, y) 穿越零点的时刻。概念清楚之后，真正落到代码里，问题就变成：给定一段数值解轨迹，如何以可接受的精度和代价，把事件时刻 t* 找出来。本节围绕四个层层递进的手段展开：时间步进扫描负责“发现”，事件函数插值负责“定位区间”，二分法与牛顿迭代负责“精化时刻”。它们不是互斥选项，而是一条典型的组合流水线。

一、时间步进扫描：用符号变化圈定事件区间

最直接的做法是在求解器输出的离散网格 t₀ < t₁ < … < tₙ 上逐点计算 gₖ = g(tₖ, yₖ)，检查相邻两点是否满足 gₖ · gₖ₊₁ ≤ 0。若成立，说明区间 [tₖ, tₖ₊₁] 内至少存在一个零点，事件被“捕获”。这一步只依赖已有的解数据，几乎零额外成本，是事件检测的第一道筛子。

但扫描有两个必须正视的陷阱。其一是漏检：若 g 在区间内先正后负再回正（两次穿越），端点同号，扫描会完全错过。其二是步长依赖：事件定位精度上限就是步长 h，粗网格下捕获的区间可能很宽。工程上的对策是：对可能高频振荡的事件函数，扫描步长应显著小于事件特征时间尺度；必要时对 g 做局部加密采样，而不是盲目缩小整个求解步长。

另一个细节是方向性。很多场景只关心上升穿越（g 由负变正）或下降穿越（g 由正变负），例如“液位首次达到上限”只取上升沿。实现时可在扫描条件中附加符号方向判断，避免同一物理事件被重复触发或误触发。

二、事件函数插值：在区间内构造可求根的近似

扫描只给出“括号区间” [tₖ, tₖ₊₁]，下一步是在区间内构造 g 的近似表达式。最简单的是线性插值：把两端点 (tₖ, gₖ) 与 (tₖ₊₁, gₖ₊₁) 连成直线，令其为零，得到

t_lin = tₖ - gₖ · (tₖ₊₁ - tₖ) / (gₖ₊₁ - gₖ)

这个公式本质上是割线法的一次迭代，计算量极小，当 g 在区间内接近线性时精度已相当可观。若求解器本身提供稠密输出（dense output）或连续扩展，还可以直接用高阶插值多项式在区间内任意点求 y(t)，再代入 g，得到比线性插值更可靠的 g(t) 近似。这就是“事件函数插值”的完整含义：不是对事件时刻插值，而是借助解的连续近似，把 g 变成区间内可反复求值的函数，为后续迭代提供基础。

三、二分法：稳健但线性收敛的保底手段

有了括号区间且 g 在端点异号，二分法就有了用武之地。每步取中点 t_m = (tₐ + t_b)/2，计算 g(t_m)：若与 g(tₐ) 同号，则令 tₐ = t_m，否则 t_b = t_m。区间长度每步减半，经过 k 步后误差上界为 (t_b - tₐ)/2ᵏ。要得到 1e-10 量级的时间精度，从 0.1 的初始区间出发约需 30 步。

二分法的价值在于无条件收敛：只要初始括号正确、g 连续，它一定收敛，不会因导数信息失真而发散。代价是收敛慢（线性收敛），且每步都要重新求值 g，若 g 的求值需要调用稠密输出，成本不可忽略。因此它常被用作“保底”和“验证”手段，而不是唯一手段。

四、牛顿迭代：快，但需要导数与好的初值

牛顿法在当前估计 tₖ 处线性化 g：tₖ₊₁ = tₖ - g(tₖ)/g′(tₖ)。若 g 光滑且初值落在零点附近，它具有二次收敛，通常 3～5 步即可达到机器精度。难点在于 g′ 的获取：解析求导往往不可行，实践中多用有限差分 g′(t) ≈ [g(t+δ) - g(t-δ)]/(2δ)，δ 取步长的平方根量级较为稳妥。

牛顿法的风险同样明确：初值离零点太远、g′ 接近零、或 g 在区间内非单调，都可能导致迭代跳出括号区间甚至发散。因此工程实现几乎总是采用混合策略：以插值结果作初值，先尝试牛顿步；若新点落在当前括号之外，或 |g| 未下降，则退回二分步。这样既保留二次收敛的速度，又继承二分法的可靠性。

五、组合策略与收敛判据

把上述环节串起来，一个可用的实现骨架大致如下（伪代码，参数为示例值）：

eps_g   = 1e-9      # 事件函数残差容差
eps_t   = 1e-10     # 时间增量容差
max_iter = 50

for k in range(n-1):
    if g[k] * g[k+1] <= 0 and direction_ok(g[k], g[k+1]):
        ta, tb = t[k], t[k+1]
        t = ta - g[k]*(tb-ta)/(g[k+1]-g[k])   # 插值初值
        for it in range(max_iter):
            gt = g_event(t)
            if abs(gt) < eps_g: break
            gp = (g_event(t+delta) - g_event(t-delta)) / (2*delta)
            tn = t - gt/gp if abs(gp) > 1e-14 else None
            if tn is None or not (ta < tn < tb) or abs(g_event(tn)) > abs(gt):
                tn = 0.5*(ta + tb)                # 退回二分
            if g_event(tn) * g_event(ta) > 0: ta = tn
            else: tb = tn
            if abs(tn - t) < eps_t: t = tn; break
            t = tn
        record_event(t, y_at(t))

几个实现要点值得强调。第一，收敛判据要双管齐下：只看 |g| < ε 在 g 平坦时可能过早停止，只看 |Δt| < εₜ 在 g 陡峭时可能精度不足，两者同时满足更稳妥。第二，括号维护不可省略：每次迭代后都要保证 ta、tb 仍夹住零点，这是混合策略可靠性的根基。第三，事件去重：若事件触发后状态发生切换，需重新初始化求解器并重置扫描起点，否则同一穿越可能被反复捕获。

六、精度与代价的权衡

四种手段的定位可以这样理解：扫描决定“能不能发现”，插值决定“初值好不好”，二分保证“一定收敛”，牛顿决定“收敛多快”。在多数工程问题中，线性插值加若干步混合迭代已能把事件时刻误差压到远小于求解器自身截断误差的水平，此时继续收紧 ε 收益有限。反过来，若事件函数本身来自实验拟合或含噪声，过度追求 1e-12 的定位精度并无意义，反而可能因 g 的微小抖动导致迭代震荡。合理的做法是先估计 g 的噪声水平与特征时间尺度，再据此设定 ε_g 与 ε_t，使数值精度与问题本身的可辨识度相匹配。