Sampling events

In this example, we will show how to use heyoka.py’s event detection system to detect when the solution of a system of ODEs satisfies certain constraints.

We will be considering the system of polynomial ODEs

\[\begin{split} \begin{cases} \frac{dx}{dt} = -y-x^2\\ \frac{dy}{dt} = 2x-y^3 \end{cases} \end{split}\]

with initial conditions

\[\begin{split} \begin{cases} x\left( 0 \right ) = 1\\ y\left( 0 \right ) = 1 \end{cases}. \end{split}\]

Let us begin by creating the symbolic variables and defining the system of ODEs:

import heyoka as hy

x, y = hy.make_vars("x", "y")
eqns = [(x, -y - x**2), (y, 2*x - y**3)]

As a first step, we want to detect when the solution crosses the \(x\) axis. The equation for this event is simply

\[ y = 0. \]

Each time the solution crosses the \(x\) axis, we want to record in a list the \(x\) coordinate at the time of crossing. Thus, we define the following callback for the event:

# The list of values of the x coordinate
# at the time of crossing.
ev_x = []

# The event's callback.
def cb(ta, time, d_sgn):
    # Determine the state of the system
    # at the time of crossing.
    ta.update_d_output(time)
    
    # Add the x coordinate at the time
    # of crossing to ev_x.
    ev_x.append(ta.d_output[0])

We can now proceed to the creation of the integrator:

ta = hy.taylor_adaptive(eqns, [1., 1.], nt_events = [hy.nt_event(
    # The event equation.
    y,
    # The callback function.
    cb)])

The event is created as a non-terminal event, as it is just a logging event that does not alter the state or the dynamics of the system.

We can then proceed to the integration of the system over a time grid:

import numpy as np
grid = np.linspace(0, 80, 2000)
out = ta.propagate_grid(grid)

Let us plot the solution and mark the \(x\) axis crossing positions:

from matplotlib.pylab import plt

fig = plt.figure(figsize=(9, 9))
ax = plt.subplot(111)
ax.set_aspect('equal', adjustable='box')

plt.plot(out[5][:, 0], out[5][:, 1])
plt.plot(ev_x, [0] * len(ev_x), 'o', alpha=.5)

plt.xlabel("x")
plt.ylabel("y")
plt.xlim((-0.75, 1))

plt.grid()

The plot shows how the event detection system indeed fired each time the solution crossed the \(x\) axis.

Let us now complicate things a bit: instead of detecting when the solution crosses the \(x\) axis, we want to detect when the solution intersects the non-linear curve

\[ y = 1.1 \sin\left( 5x \right). \]

The only thing we need to change in our setup is the event equation, which will now read

\[ y - 1.1 \sin\left( 5x \right) = 0 \]

rather than simply \(y=0\). Additionally, we will also need to store explicitly the \(y\) coordinate of the crossing points. Let us take a look at the code:

fig = plt.figure(figsize=(9, 9))

# Lists of x/y coordinates for the crossing points.
ev_x = []
ev_y = []

# The (updated) callback function.
def cb(ta, time, d_sgn):
    ta.update_d_output(time)
    ev_x.append(ta.d_output[0])
    ev_y.append(ta.d_output[1])

# Reset the integrator with the new event.
ta = hy.taylor_adaptive(eqns, [1., 1.], nt_events = [hy.nt_event(y - 1.1*hy.sin(5*x), cb)])

# Integrate over a time grid.
out = ta.propagate_grid(grid)

# Plot the results.
ax = plt.subplot(111)
ax.set_aspect('equal', adjustable='box')

plt.plot(out[5][:, 0], out[5][:, 1])
plt.plot(ev_x, ev_y, 'o', alpha = .5)

# Plot also the curve.
xrng = np.linspace(-0.75, 1, 2000)
plt.plot(xrng, 1.1*np.sin(5*xrng), 'k--')

plt.xlabel("x")
plt.ylabel("y")
plt.xlim((-0.75, 1))

plt.grid()

The plot shows again how the event detection system was able to correctly detect the intersections of the solution with the curve.