Example: Bouncing Ball

A bouncing ball is a simple hybrid system with continuous dynamics between impacts and a discrete state update at each impact.

The continuous model is:

\[\frac{d^2x}{dt^2} = -g\]

where \(x\) is height and \(g\) is gravitational acceleration. Introducing velocity \(v = dx/dt\), we write:

\[\begin{split}\begin{aligned} \frac{dx}{dt} &= v \\\\ \frac{dv}{dt} &= -g \end{aligned}\end{split}\]

with initial conditions \(x(0)=h\), \(v(0)=0\).

At each ground contact, we would ideally apply a discrete reset:

\[v^+ = -e v^-\]

where \(e\) is the coefficient of restitution.

In DiffSL we define an event (root) function stop_i { x } so integration halts when \(x = 0\). A reset_i block applies the discrete velocity reset \(v^+ = -e v^-\) at each impact. The example calls solve_dense(...) which automatically handles repeated bounce events throughout the integration, returning the full trajectory up to final_time.

import numpy as np
import matplotlib.pyplot as plt
import pydiffsol as ds


def solve():
    final_time = 10.0

    ode = ds.Ode(
        """
        in_i {
            g = 9.81,
            h = 10.0,
            e = 0.8,
        }
        u_i {
            x = h,
            v = 0.0,
        }
        F_i {
            v,
            -g,
        }
        stop {
            x,          // stop when height x reaches zero
        }
        reset_i {
            1e-12,
            -e * v,
        }
        """,
        matrix_type=ds.nalgebra_dense,
        ode_solver=ds.tsit45,
        linear_solver=ds.lu,
    )

    params = np.array([9.81, 10.0, 0.8])
    t_eval = np.linspace(0.0, final_time, 200)
    solution = ode.solve_dense(params, t_eval)

    ts = solution.ts
    x, v = solution.ys

    fig, ax = plt.subplots()
    ax.plot(ts, x, label="x (height)")
    ax.plot(ts, v, label="v (velocity)")
    ax.set_xlabel("t")
    ax.set_ylabel("state")
    ax.legend()
    fig.savefig("docs/images/bouncing_ball.svg")


if __name__ == "__main__":
    solve()