检查solve_ivp中的事件

Question

在scipy.solve_ivp()中，有什么快速的方法来检查哪个事件被触发了？
为了让大家明白，我正在模拟一个双摆。我想检查第一个杆子或第二个杆子是否翻转了。我为每个杆子创建了两个事件函数，但一旦事件被触发，有没有办法知道是哪个事件被触发了呢？谢谢 :)

Answer 1

如果你给 solve_ivp 传递多个事件，它返回的 sol.t_events 和 sol.y_events 会通过把不同的事件放在不同的数组里来告诉你发生了什么事件。

这里有一个例子，我是从 SciPy 文档 改编过来的。我在模拟一门炮的发射，我想知道炮弹什么时候达到最高点，以及什么时候落到地面。

from scipy.integrate import solve_ivp


def upward_cannon(t, y):
    return [y[1], -0.5]

def hit_ground(t, y):
    return y[0]

def peak(t, y):
    return y[1]


hit_ground.terminal = True
hit_ground.direction = -1


events = [peak, hit_ground]
sol = solve_ivp(upward_cannon, [0, 101], [0, 10], events=events)
print(sol.t_events)

输出结果：

[array([20.]), array([40.])]

在这里，events 的第一个元素是 peak()。因此，sol.t_events 的第一个元素是一个数组，里面存放的是这个事件发生的时间点。这个值是 array([20.])，所以这个事件只发生了一次，发生在 t = 20 的时候。第二个数组对应的是第二组事件。

sol.y_events 的返回值也是类似的。

Answer 2

哦，我看到已经有答案了。既然已经写好了，那我就发出来吧：

下面的代码中有两个终止事件，它表示一个简单的谐振子。

import numpy as np
from scipy import integrate

def dy(t, y):
    return [y[1], -y[0]]

def event1(t, y):
    return y[0]
event1.terminal = True

def event2(t, y):
    return y[1]
event2.terminal = True

tspan = [0, 10]
y0 = [0.1, -0.1]
sol = integrate.solve_ivp(dy, tspan, y0, events=[event1, event2])

根据solve_ivp的文档，sol.t_events是一个数组列表，里面记录了事件发生的时间。列表中的第一个数组对应第一个事件；第二个数组对应第二个事件。如果有些事件不是终止事件，可以创建一个列表，记录每个数组中最后一个事件的时间：

last_events = [(t_event[-1] if len(t_event) else tspan[0]) for t_event in sol.t_events]

发生的事件是这些数组的argmax。

np.argmax(last_events)  # 0 

在这个例子中，第一个事件（事件0）发生了，这意味着“位置”y[0]已经穿过了零点。你可以把初始的“速度”改成y0[1] = -y0[1]，这样你会看到argmax变成了1，这就意味着第二个事件（事件1）发生了，也就是说速度已经穿过了零点。