如何强制solve_ivp仅允许积分变量输出非负解？

核心约束要求 #

开展表面冰层随时间演化的冻结过程建模时，冻结/融化速率取值可正可负，需满足物理约束：当表面已冻结冰量为0时，即使冻结速率为负（即融化方向），冰量取值不得低于0，负值无物理意义。

原有写法失效原因 #

原代码在导数函数df中通过if判断直接将小于0的x[0][0]置0的操作无法生效，原因是：solve_ivp调用导数函数时传入的x是求解器步进过程中的临时状态副本，函数内部对x的修改不会回写到求解器的状态更新链路，求解器仍会按照计算得到的负导数继续推进，最终还是会输出负的冰量结果。

可行实现方案 #

方案1：使用solve_ivp内置非负约束参数（最简便） #

solve_ivp原生支持指定状态量的非负约束，无需手动修改导数逻辑，仅需在调用求解器时传入non_neg参数标记需要约束非负的状态即可，LSODA求解器原生兼容该参数。
修正后的代码示例：

def df(t,x,*system_constants):
    dxdt = np.zeros_like(x)
    # 计算原始冻结/融化速率
    dxdt[0][0] = x[0][0] * a
    return dxdt

# 初始条件与时间范围配置
t0 = 0
tf = 10 
t_span = (t0, tf)

# 配置非负约束：标记冰量对应状态位不可为负
non_neg_mask = np.zeros_like(x0, dtype=bool)
non_neg_mask[0, 0] = True

# 求解传入非负约束参数
sol = solve_ivp(
    df, 
    t_span, 
    x0, 
    args=system_constants, 
    method='LSODA', 
    dense_output=True,
    non_neg=non_neg_mask
)

传入该参数后，求解器会在步进过程中自动检测被标记状态的取值，触碰到0边界时自动调整步长、修正状态更新逻辑，保证冰量不会降到0以下。

方案2：导数层边界速率截断（适合自定义边界逻辑场景） #

如果需要自定义边界触发后的行为，不要修改输入状态x本身，而是在冰量触达0边界时直接截断对应导数，避免出现驱动冰量继续下降的负速率：

def df(t,x,*system_constants):
    dxdt = np.zeros_like(x)
    # 计算原始速率
    raw_rate = x[0][0] * a
    # 边界修正：冰量接近0且速率为融化方向时，速率置0
    if x[0][0] <= 1e-12 and raw_rate < 0:
        dxdt[0][0] = 0.0
    else:
        dxdt[0][0] = raw_rate
    return dxdt

判断时加入1e-12的浮点容差，避免求解器数值计算误差导致边界判断失效。该方案的物理逻辑是：冰量耗尽后无冰可融，因此融化方向的速率直接截断为0，自然不会产生负冰量。

注意事项 #

• 禁止在导数函数内部直接修改输入的状态变量x，该操作不会被求解器识别，属于无效操作
• 两种方案二选一即可，不要同时使用内置非负约束和手动速率截断，避免双重修正引入求解误差
• 正式使用前建议做简单验证：设置恒定负融化速率，分别取初始冰量为0、初始冰量为正值两组初始条件，确认冰量为0时不会继续下降、冰量降到0后能稳定保持在0值。

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