Checking Event in solve_ivp

Question

What is a fast way of checking which event has been triggered in scipy.solve_ivp()? For context, I am simulating a double pendulum. I want to check if either the first rod or the second rod has been flipped. I have created two event functions for each, but once the event is triggered, is there any way I can retrieve which event has been triggered? Thanks :)

Answer 1

Oops, I see there's an answer. Might as well post, since it's already written:

There are two terminal events in the following code, which represents a simple harmonic oscilator.

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])

According to the solve_ivp documentation, sol.t_events is a list of arrays of times at which an event occured. The first array in the list corresponds with the first event; the second array corresponds with the second event. In case some of the events are non-terminal, create a list with time of the last event in each array:

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

The event which occured is the argmax of these.

np.argmax(last_events)  # 0 

In this case, the first event (event 0) occured, which means that the "position" y[0] has crossed zero. You can change the initial "velocity" to y0[1] = -y0[1] to see that the argmax becomes 1, so the second event (event 1) occurs, which means that the velocity has crossed zero.

Answer 2

If you pass multiple events to solve_ivp, the return value sol.t_events and sol.y_events indicates which kind of event happened by putting different events in different arrays.

Here's an example, which I have adapted from the SciPy docs. I am simulating a cannon being fired, and I want to know both when the cannonball reaches its peak, and when the cannonball hits the ground.

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)

Output:

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

Here, the first element in events is peak(). Because of that, the first element of sol.t_events is an array of t values where that event occurred. That value is array([20.]), so it happened only once, at t = 20. The second array corresponds to the second set of events.

The sol.y_events return value is similar.