How to do ODE events

15.1 Introduction

Solving differential equations, we are often interesting in knowing if certains events occur and if so at what time. For instance, if you solve a differential equation for a ball falling in a gravitational field, you might want to know the time for which the ball hits the ground, and stop the integration of the differential equation at that time.

Using solve_ivp it is possible to provide extra functions, event-functions, that monitor mathematically if any events occur. Each function must take the same input as the dydt-function (the first argument of solve_ivp), and must return a float which is a continuous value being a function of the independent and dependent variables in the problem. Given one or several such event-functions, solve_ivp will monitor if there are any sign-changes in the floats of the event-functions, while it solves the differential equations. If so, the precise times at which that happens will be evaluated and returned in the end. If the event-function whose float is changing sign has an attribute, .terminal set to True, then solve_ivp will further stop at the time of the event.

15.2 Collect values

As a first example consider that we are solving:

$\frac{dy}{dt}=0.9 y^2\cos t,\quad (t,y)=(0,1)$

and are interested in knowing at which times, $t$ , the variable $y$ become $5$ . To do this we have to introduce an event function that changes sign continuously whenever $y$ becomes $5$ . An obvious choice would be:

$f(y) = y-5,$

but other functions would be equally good, e.g.:

$f_\mathrm{alt-1}(y) = 5-y,$

$f_\mathrm{alt-2}(y) = \exp(y-5) -1.$

Now using the simplest choice, $f$ , it looks like this:

%matplotlib inline
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
plt.rc('font', size=16)

tinit = 0
tfinal = 10
trange = [tinit,tfinal]
t_eval = np.linspace(tinit,tfinal,100)
yinit = np.array([1]) # the event function requires it to be np.array

dydt = lambda t, y: 0.9 * y**2 * np.cos(t)

# define the event function with the
# independent and dependent variable as input
f = lambda t, y: y - 5

# supply it in a list as value for keyword "events":
mysol = solve_ivp(dydt, trange, yinit, t_eval=t_eval, max_step=1e-2,
                  events=[f])
ts = mysol.t
ys = mysol.y[0]

# now extract the times at which the events occurred
t_events = mysol.t_events[0] # [0] because we want output from the
                             # first (and only) event function
print(t_events)

fig,ax = plt.subplots(1,1)
ax.plot(ts,ys)
ax.grid()
ax.set_ylim([0,11])
ax.set_xlabel('$t$')
ax.set_ylabel('$y$')

ax.scatter(t_events,5 * np.ones(t_events.shape))

fig.tight_layout()
fig.savefig('ode_events_1.png')

[1.09491408 2.04667858 7.37809938 8.32986388]

Figure 1

15.3 Two event functions

Here we add one more event function that monitors whenever the condition:

$y = t$

if fulfilled.

%matplotlib inline
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
plt.rc('font', size=16)

tinit = 0
tfinal = 10
trange = [tinit,tfinal]
t_eval = np.linspace(tinit,tfinal,100)
yinit = np.array([1])

dydt = lambda t, y: 0.9 * y**2 * np.cos(t)

# define the event functions 
f = lambda t, y: y - 5
g = lambda t, y: y - t

# supply it in a list as value for keyword "events":
mysol = solve_ivp(dydt, trange, yinit, t_eval=t_eval, max_step=1e-2,
                  events=[f,g]) # provide two event functions
ts = mysol.t
ys = mysol.y[0]

t_events_f = mysol.t_events[0] # times for events with f
t_events_g = mysol.t_events[1] # times for events with g

fig,ax = plt.subplots(1,1)
ax.plot(ts,ys)
ax.grid()
ax.set_ylim([0,11])
ax.set_xlabel('$t$')
ax.set_ylabel('$y$')

print(t_events)
ax.scatter(t_events_f,5 * np.ones(t_events_f.shape))
ax.plot(t_events_g,t_events_g,marker='s',linestyle='dashed')

fig.tight_layout()
fig.savefig('ode_events_2.png')

Figure 2

15.4 Stop integration

When event functions are defined using def, the resulting function may be given any of the two attributes .direction and .terminal.

If the .direction attribute is set to some positive value, it means that solve_ivp should only consider an event to occur if the event function changes its value in ascending direction (from negative to positive). Likewise, if the .direction attribute is set to some negative value, it means that an event has only occured whenever the event function changes in descending direction (from positive to negative).

If the .terminal attribute is set, solve_ivp should stop its integration whenever the first event occurs according to this event function.

Here is an example:

%matplotlib inline
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
plt.rc('font', size=16)

tinit = 0
tfinal = 10
trange = [tinit,tfinal]
yinit = np.array([1])

dydt = lambda t, y: 0.9 * y**2 * np.cos(t)

# define an event function 
def f(t,y):
    return y - 5

f.direction = -1    # only event if f changes sign while decreasing
f.terminal = True   # solve_ivp should stop

mysol = solve_ivp(dydt, trange, yinit, max_step=1e-2,
                  events=[f])
ts = mysol.t
ys = mysol.y[0]

t_events = mysol.t_events[0] # times for events with f

fig,ax = plt.subplots(1,1)
ax.plot(ts,ys)
ax.grid()
ax.set_ylim([0,11])
ax.set_xlabel('$t$')
ax.set_ylabel('$y$')

print(t_events)
ax.scatter(t_events,5 * np.ones(t_events.shape))

fig.tight_layout()
fig.savefig('ode_events_3.png')

And the differential equation is now only solved up to the point that the event function,

$f(y) = y - 5$

is changing from positive to negative for the first time, as evidenced by the resulting figure:

Figure 3