Event Detection with CVODE¶

SUNDALS CVODE includes a rootfinding feature: given event functions $g_i(t, y)$, the solver detects zero crossings to high accuracy without requiring the user to manually bracket or bisect.

Pass events as the 8th argument in expression mode:

np_cvode(rhs, vars, y0, tspan, rtol, atol, method, [g1, g2, ...])

When events are present, np_cvode returns [trajectory, events_list] where each event is [t_event, [y_values], [indices]].

In [11]:
load("numerics")$
load("numerics-sundials")$
load("ax-plots")$

1. Bouncing Ball¶

Free-fall with gravity: $\dot{h} = v$, $\dot{v} = -g$.

Event: $h = 0$ (ball hits ground). We detect the impact, then could use the event time and velocity to restart with a coefficient of restitution.

Drop from $h = 10\,$m, $v_0 = 0$. Expected impact at $t = \sqrt{2h/g} \approx 1.43\,$s.

In [12]:
g : 9.81$
tspan1 : np_linspace(0.0, 3.0, 300)$

[traj1, events1] : np_cvode(
    [v, -g], [t, h, v], [10.0, 0.0], tspan1,
    1e-12, 1e-12, adams, [h])$

printf(true, "Number of events: ~d~%", length(events1))$
ev1 : first(events1)$
printf(true, "Impact time: ~,6f s  (exact: ~,6f)~%",
       first(ev1), sqrt(2*10.0/g))$
printf(true, "Impact velocity: ~,4f m/s~%", second(ev1)[2])$

Number of events: 1
o114Impact time: 1.427843 s  (exact: 1.427843)
o114Impact velocity: -14.0071 m/s
o114$$\mathbf{false}$$
false
In [13]:
t1 : np_slice(traj1, all, 0)$
h1 : np_slice(traj1, all, 1)$

ax_draw2d(
  color="#2563eb", line_width=2, name="h(t)",
  lines(t1, h1),
  color="#dc2626", marker_size=8, name="impact",
  points([first(ev1)], [second(ev1)[1]]),
  title="Bouncing Ball (single bounce)",
  xlabel="t (s)", ylabel="h (m)",
  grid=true, showlegend=true
)$
No description has been provided for this image

2. Bouncing Ball — Multiple Bounces¶

To simulate repeated bounces, we restart the solver after each impact, reversing the velocity with a coefficient of restitution $e = 0.8$.

In [14]:
/* Multi-bounce simulation via restart */
cor : 0.8$  /* coefficient of restitution */
t_end : 8.0$
dt : 0.01$

h_all : []$  v_all : []$  t_all : []$
event_times : []$
t_start : 0.0$  h0 : 10.0$  v0 : 0.0$

for bounce : 1 thru 20 do (
  if t_start >= t_end then return(done),
  ts : makelist(t_start + i*dt, i, 0, ceiling((t_end - t_start)/dt)),
  ts : sublist(ts, lambda([x], x <= t_end)),
  if length(ts) < 2 then return(done),
  sol : np_cvode([bv, -g], [bt, bh, bv],
                 [h0, v0], ts, 1e-12, 1e-12, adams, [bh]),
  /* np_cvode returns [traj, events] when events param is given */
  traj : first(sol),
  evs : second(sol),
  /* Collect trajectory */
  n : first(np_shape(traj)),
  for i : 0 thru n-1 do (
    t_all : endcons(np_ref(traj, i, 0), t_all),
    h_all : endcons(np_ref(traj, i, 1), h_all),
    v_all : endcons(np_ref(traj, i, 2), v_all)
  ),
  if length(evs) = 0 then return(done),
  ev : first(evs),
  event_times : endcons(first(ev), event_times),
  t_start : first(ev) + 1e-10,
  h0 : 0.0,
  v0 : -cor * second(ev)[2]
)$

printf(true, "~d bounces detected~%", length(event_times))$

4 bounces detected
o130$$\mathbf{false}$$
false
In [15]:
ax_draw2d(
  color="#2563eb", line_width=1.5,
  lines(t_all, h_all),
  title="Bouncing Ball with Restitution (e=0.8)",
  xlabel="t (s)", ylabel="h (m)",
  grid=true
)$
No description has been provided for this image

3. Projectile with Drag — Target Altitude¶

A projectile launched upward with quadratic air drag: $$\dot{h} = v, \quad \dot{v} = -g - \frac{\rho\,C_d\,A}{2m}\,v\,|v|$$

Events: detect when the projectile crosses $h = 50\,$m (ascending and descending).

In [16]:
/* Physical parameters */
mass : 1.0$    /* kg */
Cd : 0.47$     /* drag coefficient (sphere) */
Acs : 0.01$    /* cross-section area m^2 */
rho : 1.225$   /* air density kg/m^3 */
drag_coeff : rho * Cd * Acs / (2 * mass)$

/* v*abs(v) for correct drag direction */
tspan3 : np_linspace(0.0, 12.0, 600)$

[traj3, events3] : np_cvode(
    [pv, -g - drag_coeff * pv * abs(pv)],
    [pt, ph, pv], [0.0, 50.0], tspan3,
    1e-10, 1e-10, adams, [ph - 50.0])$

printf(true, "Events at h=50m:~%")$
for ev in events3 do
  printf(true, "  t = ~,4f s, h = ~,4f m, v = ~,4f m/s~%",
         first(ev), second(ev)[1], second(ev)[2])$

Events at h=50m:
o149  t = 1.2329 s, h = 50.0000 m, v = 31.9721 m/s
o149  t = 7.3296 s, h = 50.0000 m, v = -28.0417 m/s
o149$$\mathbf{done}$$
false
In [17]:
t3 : np_slice(traj3, all, 0)$
h3 : np_slice(traj3, all, 1)$

ev_t3 : map(first, events3)$
ev_h3 : map(lambda([e], second(e)[1]), events3)$

ax_draw2d(
  color="#2563eb", line_width=2, name="h(t)",
  lines(t3, h3),
  color="#dc2626", marker_size=8, name="h=50m crossings",
  points(ev_t3, ev_h3),
  color="#9ca3af", dash="dash",
  explicit(50, t, 0, 12),
  title="Projectile with Drag",
  xlabel="t (s)", ylabel="h (m)",
  grid=true, showlegend=true
)$
No description has been provided for this image

4. Damped Pendulum — Zero-Crossing Detection¶

A damped pendulum: $$\ddot{\theta} + 2\zeta\omega_0\,\dot{\theta} + \omega_0^2\sin\theta = 0$$

Event: $\theta = 0$ (each time the pendulum passes through the vertical). We detect all zero-crossings and measure the decaying amplitude.

In [18]:
omega0 : 2*float(%pi)$  /* natural frequency: 1 Hz */
zeta : 0.05$             /* light damping */

tspan4 : np_linspace(0.0, 10.0, 2000)$

[traj4, events4] : np_cvode(
    [w, -2*zeta*omega0*w - omega0^2*sin(th)],
    [t, th, w], [float(%pi/4), 0.0], tspan4,
    1e-12, 1e-12, adams, [th])$

printf(true, "Zero crossings of theta:~%")$
for ev in events4 do
  printf(true, "  t = ~,4f s, omega = ~,4f rad/s~%",
         first(ev), second(ev)[2])$

Zero crossings of theta:
o169  t = 0.2688 s, omega = -4.4271 rad/s
o169  t = 0.7836 s, omega = 3.7745 rad/s
o169  t = 1.2944 s, omega = -3.2202 rad/s
o169  t = 1.8023 s, omega = 2.7484 rad/s
o169  t = 2.3082 s, omega = -2.3465 rad/s
o169  t = 2.8126 s, omega = 2.0038 rad/s
o169  t = 3.3160 s, omega = -1.7115 rad/s
o169  t = 3.8186 s, omega = 1.4619 rad/s
o169  t = 4.3207 s, omega = -1.2489 rad/s
o169  t = 4.8224 s, omega = 1.0669 rad/s
o169  t = 5.3238 s, omega = -0.9116 rad/s
o169  t = 5.8250 s, omega = 0.7788 rad/s
o169  t = 6.3261 s, omega = -0.6654 rad/s
o169  t = 6.8270 s, omega = 0.5686 rad/s
o169  t = 7.3278 s, omega = -0.4858 rad/s
o169  t = 7.8286 s, omega = 0.4151 rad/s
o169  t = 8.3294 s, omega = -0.3547 rad/s
o169  t = 8.8301 s, omega = 0.3031 rad/s
o169  t = 9.3308 s, omega = -0.2590 rad/s
o169  t = 9.8314 s, omega = 0.2213 rad/s
o169$$\mathbf{done}$$
false
In [19]:
t4 : np_slice(traj4, all, 0)$
th4 : np_slice(traj4, all, 1)$

ev_t4 : map(first, events4)$

ax_draw2d(
  color="#2563eb", line_width=1.5, name="theta(t)",
  lines(t4, th4),
  color="#dc2626", marker_size=6, name="zero crossings",
  points(ev_t4, makelist(0, i, 1, length(ev_t4))),
  color="#9ca3af", dash="dot",
  explicit(float(%pi/4) * exp(-zeta*omega0*t), t, 0, 10),
  title="Damped Pendulum",
  xlabel="t (s)", ylabel="theta (rad)",
  grid=true, showlegend=true
)$
No description has been provided for this image

5. Multiple Simultaneous Events¶

Monitor two event functions at once: a linear ramp $\dot{y} = 1$ crossing thresholds at $y = 0.3$ and $y = 0.7$.

In [20]:
[traj5, events5] : np_cvode(
    [1], [t, y], [0.0], [0.0, 0.5, 1.0],
    1e-12, 1e-12, adams, [y - 0.3, y - 0.7])$

for ev in events5 do
  printf(true, "t = ~,6f, y = ~,6f, event indices = ~a~%",
         first(ev), second(ev)[1], third(ev))$

t = 0.300000, y = 0.300000, event indices = [0]
o185t = 0.700000, y = 0.700000, event indices = [1]
o185$$\mathbf{done}$$
false