能否让scipy.integrate.solve_ivp拒绝步长以避免无效状态下计算RHS？

问题场景 #

当使用solve_ivp求解带终端事件的微分方程时，求解器可能会跨大步越过设定的事件阈值，导致微分方程函数收到无效参数（如示例中y<0）。示例代码如下：

import numpy as np
from scipy.integrate import solve_ivp

def funct(t, y):
    return -np.sqrt(y)

def event(t, y):
    return y[0]-0.1

if __name__ == '__main__':
    event.terminal = True
    sol = solve_ivp(funct, t_span=[0, 3], y0=[1], events=event)

示例中仅触发警告，但实际场景中若调用第三方库，无效参数会直接导致求解失败。

解决方案 #

1. 在微分方程函数中加入安全阈值检查 #

直接在函数内判断参数合法性，当接近阈值时主动抛出异常，迫使求解器缩小步长：

def funct(t, y):
    safe_threshold = 0.05  # 设为比事件阈值稍低的安全值
    if y[0] < safe_threshold:
        raise ValueError("y值低于安全阈值，步长过大")
    return -np.sqrt(y)

同时配合max_step参数限制最大步长，降低越过阈值的概率：

sol = solve_ivp(funct, t_span=[0, 3], y0=[1], events=event, max_step=0.01)

2. 优化事件函数的方向检测 #

为事件函数添加方向标记，让求解器更精准识别事件触发的方向，减少越过阈值的情况：

event.direction = -1  # 表示当y从高到低穿过阈值时触发事件

3. 动态调整步长的回调函数 #

利用on_step回调，在每一步计算后检查参数，接近阈值时强制缩小下一步的步长：

def on_step(t, y):
    if y[0] < 0.2:  # 提前预警，阈值设为事件阈值的1.2~2倍
        sol.step.max_step = 0.001  # 动态调小最大步长
    return True

# 调用时传入回调
sol = solve_ivp(funct, t_span=[0, 3], y0=[1], events=event, on_step=on_step)

4. 切换到事件检测更精准的求解器 #

隐式求解器（如Radau、BDF）对终端事件的处理精度通常优于显式求解器（如RK45），可尝试切换求解器：

sol = solve_ivp(funct, t_span=[0, 3], y0=[1], events=event, method='Radau')

内容的提问来源于stack exchange，提问作者mauro