Event Detection Example — Bouncing Ball

This example demonstrates zero-crossing (event) detection using differentialequations.ivp_with_events(). A ball is dropped from 10 m and bounces under gravity with a coefficient of restitution of 0.85. Each floor impact (height = 0) is detected as a terminal event and the integration is restarted with the reflected velocity for up to 5 bounces.

\[y'' = -g, \quad g = 9.81\;\mathrm{m/s^2}\]

State vector: \([h, v]\) (height, velocity).

Example Code

"""IMSL event detection example: Bouncing ball simulation.

Demonstrates the use of ivp_with_events() for zero-crossing (event) detection.
Simulates a bouncing ball under gravity with coefficient of restitution 0.85:

    y'' = -g   (g = 9.81 m/s²)

State: y = [height, velocity].  Event: height == 0 (floor impact).
The simulation runs for 5 bounces by restarting integration after each impact.

Outputs:
- Impact times and velocities printed to stdout
- SVG plot saved to test_output/example_imsl_events.svg
"""

from __future__ import annotations

from pathlib import Path
from typing import Dict, List

import matplotlib.pyplot as plt
import numpy as np

from differentialequations import ivp_with_events


_G = 9.81
_COR = 0.85
_MAX_BOUNCES = 5


def _gravity_rhs(t: float, y: np.ndarray) -> list:
    """Gravity ODE right-hand side.

    Args:
        t (float): Time (unused, autonomous).
        y (np.ndarray): State [height, velocity].

    Returns:
        list: [velocity, acceleration] = [y[1], -g].
    """
    return [y[1], -_G]


def _floor_event(t: float, y: np.ndarray) -> float:
    """Event function: returns height (zero-crossing = floor impact).

    Args:
        t (float): Time.
        y (np.ndarray): State [height, velocity].

    Returns:
        float: Height above floor (event at y[0] == 0).
    """
    return y[0]


_floor_event.terminal = True
_floor_event.direction = -1


def run_demo_imsl_events() -> Dict[str, object]:
    """Run IMSL event detection example: Bouncing ball with 5 bounces.

    Integrates a bouncing ball under gravity using ivp_with_events() to detect
    floor impacts (height == 0) and restarts integration after each bounce,
    applying coefficient of restitution (COR = 0.85) to reverse velocity.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``impact_times`` (list),
            ``impact_velocities`` (list), and ``plot_path`` (str).
    """
    y0 = [10.0, 0.0]
    t_start = 0.0
    t_max_segment = 30.0

    all_t: List[np.ndarray] = []
    all_y: List[np.ndarray] = []
    impact_times: List[float] = []
    impact_velocities: List[float] = []

    y_current = list(y0)
    t_current = t_start

    for bounce in range(_MAX_BOUNCES):
        t_end = t_current + t_max_segment
        result = ivp_with_events(
            _gravity_rhs,
            (t_current, t_end),
            y_current,
            events=[_floor_event],
            method="RK45",
        )
        all_t.append(result.t)
        all_y.append(result.y)

        if result.t.size >= 2:
            t_impact = result.t[-1]
            v_impact = result.y[1, -1]
            impact_times.append(t_impact)
            impact_velocities.append(v_impact)

        v_after = -_COR * result.y[1, -1]
        if abs(v_after) < 0.01:
            break
        y_current = [0.0, v_after]
        t_current = result.t[-1]

    t_full = np.concatenate(all_t)
    y_full = np.concatenate(all_y, axis=1)

    print("\nIMSL Event Detection Example: Bouncing Ball (COR = 0.85)")
    print("=" * 58)
    print(f"  Initial height   : {y0[0]:.2f} m")
    print(f"  Initial velocity : {y0[1]:.2f} m/s")
    print(f"  Gravity          : {_G} m/s²")
    print(f"  Restitution      : {_COR}")
    print()
    print(f"  {'Bounce':<8} {'Impact time (s)':>18} {'Impact velocity (m/s)':>22}")
    print("  " + "-" * 50)
    for i, (ti, vi) in enumerate(zip(impact_times, impact_velocities), start=1):
        print(f"  {i:<8} {ti:>18.4f} {vi:>22.4f}")
    print("=" * 58)

    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_events.svg"

    fig, axes = plt.subplots(1, 2, figsize=(13, 5))

    axes[0].plot(t_full, y_full[0], color="#0369a1", linewidth=1.8, label="Height")
    for ti in impact_times:
        axes[0].axvline(ti, color="#e11d48", linewidth=0.8, linestyle="--", alpha=0.6)
    axes[0].set_xlabel("Time (s)")
    axes[0].set_ylabel("Height (m)")
    axes[0].set_title("Bouncing Ball: Height vs Time")
    axes[0].set_ylim(bottom=-0.1)
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)

    axes[1].bar(
        range(1, len(impact_velocities) + 1),
        [-v for v in impact_velocities],
        color="#0891b2",
        edgecolor="#0369a1",
    )
    axes[1].set_xlabel("Bounce number")
    axes[1].set_ylabel("|Impact velocity| (m/s)")
    axes[1].set_title("Impact Speed at Each Bounce")
    axes[1].set_xticks(range(1, len(impact_velocities) + 1))
    axes[1].grid(True, axis="y", alpha=0.3)

    fig.suptitle("IMSL ivp_with_events: Bouncing Ball (COR=0.85)", fontsize=13, fontweight="bold")
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    return {
        "impact_times": impact_times,
        "impact_velocities": impact_velocities,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_imsl_events()

Plot Output

Bouncing ball height vs time and impact speed bar chart

Left: height versus time with vertical dashed lines at each impact. Right: impact speed at each successive bounce — exponential decay at rate \(\text{COR}^n = 0.85^n\).

Console Output

IMSL Event Detection Example: Bouncing Ball (COR = 0.85)
==========================================================
  Initial height   : 10.00 m
  Initial velocity : 0.00 m/s
  Gravity          : 9.81 m/s┬▓
  Restitution      : 0.85

  Bounce      Impact time (s)  Impact velocity (m/s)
  --------------------------------------------------
  1                    1.4278               -14.0071
  2                    3.8552               -11.9061
  3                    5.9184               -10.1202
  4                    7.6722                -8.6021
  5                    9.1628                -7.3118
==========================================================