n-th crossing with event detection in scipy.integrate.solve_ivp

Question

I would like to know if there is any way to stop the ODE integration after a specific number of detected events. I know that scipy.integrate.solve_ivp allows making an event terminate or not, but I want to know if it’s possible to make it terminate at the n-th event, where n is a integer number bigger than one.

For example, is it possible to compute only the first three interactions with an event and then stop the integration? If this is not possible, is there an alternative library to do it?

I read the official documentation of scipy related with the integration of ODE but I haven't found anything.

Answer 1

Here is an example hack. As suggested in the comments, wrapped_event returns an arbitrary value before the required number of events has been reached, then it starts returning the true event function. (There is a little extra complexity due to sign issues when detecting event crossings.) All you should need to do to use it with your own code is to wrap your event function with wrap_event, that pass the result into solve_ivp. Respecting the terminal attribute of the event is not requested, but easy to accommodate.

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

def wrap_event(event, n_events, t0, y0):
    def wrapped_event(t, y, n_events=n_events):
        new_event = event(t, y)
        if wrapped_event.count < n_events:
            new_sign = np.sign(new_event)
            if new_sign == -wrapped_event.sign:
                wrapped_event.count += 1
                wrapped_event.sign = new_sign
            return -1
        return new_event * (-1)**(n_events)
    wrapped_event.count = 1
    wrapped_event.sign = np.sign(event(t0, y0))
    wrapped_event.terminal = event.terminal
    return wrapped_event

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

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

t0 = 0
tf = 100
y0 = [1, 0]
n_events = 3

sol = integrate.solve_ivp(f, (t0, tf), y0, dense_output=True, 
                          events=wrap_event(event, n_events, t0, y0))
t = np.linspace(sol.t[0], sol.t[-1], 300)
x = sol.sol(t)
plt.plot(t, x[0])