在自定义 predicate 上终止Scipy solve_ivp

9gm1akwq  于 2022-11-10  发布在  其他
关注(0)|答案(1)|浏览(131)

我有一个常微分方程dy/dt = f(y,t),其中y是一个N维向量,我想用scipy.integrate.solve_ivp函数来求解。
然而,在一些情况下,如果某个 predicate g(y,t)的计算结果为True,我希望停止积分。这里的用例是,我希望y的值在积分持续时间t_end结束之前收敛到某个常量值y0。我对这个常数值y0感兴趣,并且希望通过一旦收敛就终止积分来保存时间。
我希望我可以创建一个数组来存储最后5个积分步骤中y的值,如果它们非常接近,则认为已经发生了收敛。
solve_ivpevent函数在我的例子中没有什么帮助:没有我希望找到的根,而且我对收敛发生时的t不感兴趣。我很惊讶,这种寻找收敛的看似“常见”的用例并不容易完成,而且我在Stackoverflow上也找不到类似的问题。
如果有人有什么想法,我很乐意听听。

sg24os4d

sg24os4d1#

这是访问solve_ivp在后台使用的积分器类的一个很好的候选对象。如果我们采用初始条件为y(0) = 100的简单函数dy/dt = -y。我们希望在模拟的1秒内解的变化小于0.1时终止函数,即abs(y(t) - y(t-1)) < 0.1。对于这个常微分方程,这种情况发生在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

相关问题