scipy.solve_ivp:在第二次发生事件时终止轨迹

scipy.solve_ivp: Terminate trajectory on second occurrence of event

我想求解 return 到初始条件所需的时间。作为玩具示例,让我们以一个简单的谐波振荡器 d/dt (x, y) = (y, -x) 为例,初始条件 (x,y) = (1,0).

我希望代码为 return T = 6.28... (= 2π),因为这是 return 到 y 为正的初始条件所花费的时间。当达到这一点时,它应该停止积分。

import numpy as np
from scipy.integrate import solve_ivp

def f(t,x):
    return [x[1], -x[0]]

def event(t,x):
    return x[0]

event.direction = 1
event.terminal = True # Try false, this gives the result but integrates beyond the event I need.

initial_condition = [0, 1]

sol = solve_ivp(f, [0, 10], initial_condition, method='RK45',dense_output=True, \
                events = event, rtol=1e-6, atol=3e-6)

print(sol.t_events)

问题是设置 event.terminal = True 会立即启动积分,因为事件在 t = 0 时得到满足。有没有什么办法可以让积分仅在事件满足两次后停止?

为了将来参考,这里是我一直在使用的解决方法:

将初始条件设置为事件前的一小部分,例如initial_condition = [1e-100, 1] 而不是 [0, 1]。它似乎不会影响轨迹,即使对于更高的 degree-of-freedom 混沌系统也是如此,但它不会在 t=0 时触发事件。

另一种可行的解决方法是在没有事件的情况下集成一段时间,然后使用结果作为进一步集成的初始条件,这次有事件。我怀疑这个答案对任何人都有用,但我发布这个只是为了完整性。

与您的解决方法类似,我将 1e-100 添加到事件而不是初始状态。这样做的好处是完全不影响轨迹,甚至在理论上也不影响。

当然,这确实意味着您未来的事件会提前一小部分发生,但如果这对您很重要,则取决于您的情况。

对于您编写的示例代码,这将使您的事件如下:

def event(t,x):
    return x[0] + 1e-100

我试过的一件事 没有 工作是在集成期间修改终止参数。我想如果我从 event.terminal = False 开始,然后在事件定义集 event.terminal = True 内,那么它将在第一个事件期间(在 t=0)继续,并在第二个事件处终止。但是,如果在集成过程中更改了参数,它似乎会忽略更改的参数。

(此问题与 相关,但我没有标记任何内容的名誉。)