在自定义谓词上终止 Scipy solve_ivp答案

【问题标题】：Terminate Scipy solve_ivp on custom predicate在自定义谓词上终止 Scipy solve_ivp
【发布时间】：2022-06-25 19:56:27
【问题描述】：

我有一个 ODE dy/dt = f(y,t)，其中 y 是一个 N 维向量，我想使用 scipy.integrate.solve_ivp 函数求解。

但是，如果某个谓词 g(y,t) 的计算结果为 True，我想停止集成。我在这里的用例是我希望y 的值在积分持续时间t_end 结束之前收敛到某个恒定值y0。我对这个常数值y0 感兴趣，并希望在收敛后终止集成以节省时间。

我希望在最后 5 个集成步骤中创建一个数组来存储 y 的值，如果它们非常接近，则认为已经发生收敛。

solve_ivp 的 event 函数在我的情况下并没有真正的帮助：没有我希望找到的根，并且当收敛发生时我对 t 不感兴趣。令我惊讶的是，这种寻找收敛的看似“常见”的用例并不容易完成，而且我在 Stackoverflow 上也找不到类似的问题。

如果有人有什么想法，我很乐意听到。

【问题讨论】：

为什么events 参数不起作用？你不能修改g 以返回False 而不是True（当你希望它停止时它实际上会返回零）

标签： numpy scipy

【解决方案1】：

这是访问 solve_ivp 在后台使用的集成器类的理想选择。如果我们采用具有初始条件y(0) = 100 的简单函数dy/dt = -y。当解在模拟的 1 秒内变化小于 0.1 时，我们希望终止函数，即abs(y(t) - y(t-1)) < 0.1。对于此 ODE，这发生在 t=-ln(0.1 / (100(e-1)) 或 t~7.45。我们可以使用 RK45 积分器 (RK45 docs) 解决这个问题，如下所示：

import numpy as np
from scipy.integrate import RK45
def fun(t, y):
    return -y
    
y0 = [100]
t0 = 0

#max_step used to ensure that we take small enough time steps
rk45 = RK45(fun, t0, y0, t_bound=1000, max_step=0.1)
t = []
y = []
while rk45.status == "running":
    t.append(rk45.t)
    y.append(rk45.y[0])
    if rk45.t > 1.0 and np.abs(np.interp(rk45.t-1, t, y) - rk45.y[0]) < 0.1:
        break
    rk45.step()

print(f"Final t: {t[-1]:.1f}")

# Because max_step=0.1, t[-11] will be 1 second behind t[-1]
print(f"Time period checked: {t[-1]-t[-11]:.1f}, delta_y: {y[-11]-y[-1]:.1f}")

产量

Final t: 7.5
Time period checked: 1.0, delta_y: 0.1

【讨论】：