Example: Bouncing Ball
Modelling a bouncing ball is a simple and intuitive example of a system with discrete events. The ball is dropped from a height \(h\) and bounces off the ground with a coefficient of restitution \(e\). When the ball hits the ground, its velocity is reversed and scaled by the coefficient of restitution, and the ball rises and then continues to fall until it hits the ground again. This process repeats until halted.
The second order ODE that describes the motion of the ball is given by:
\[ \frac{d^2x}{dt^2} = -g \]
where \(x\) is the position of the ball, \(t\) is time, and \(g\) is the acceleration due to gravity. We can rewrite this as a system of two first order ODEs by introducing a new variable for the velocity of the ball:
\[ \begin{align*} \frac{dx}{dt} &= v \\ \frac{dv}{dt} &= -g \end{align*} \]
where \(v = \frac{dx}{dt}\) is the velocity of the ball. This is a system of two first order ODEs, which can be written in vector form as:
\[ \frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y}, t) \]
where
\[ \mathbf{y} = \begin{bmatrix} x \\ v \end{bmatrix} \]
and
\[ \mathbf{f}(\mathbf{y}, t) = \begin{bmatrix} v \\ -g \end{bmatrix} \]
The initial conditions for the ball, including the height from which it is dropped and its initial velocity, are given by:
\[ \mathbf{y}(0) = \begin{bmatrix} h \\ 0 \end{bmatrix} \]
When the ball hits the ground, we need to update the velocity of the ball according to the coefficient of restitution, which is the ratio of the velocity after the bounce to the velocity before the bounce. The velocity after the bounce \(v’\) is given by:
\[ v’ = -e v \]
where \(e\) is the coefficient of restitution. We will first implement this using Diffsol’s procedural event-handling API, then rewrite the same event declaratively using the DiffSL stop and reset tensors.
Procedural approach
To implement the bounce in Rust, we need to detect when the ball hits the ground. We can do this by using Diffsol’s event handling feature, which allows us to specify a function that is equal to zero when the event occurs, i.e. when the ball hits the ground. This function \(g(\mathbf{y}, t)\) is called an event or root function, and for our bouncing ball problem, it is given by:
\[ g(\mathbf{y}, t) = x \]
where \(x\) is the position of the ball. When the ball hits the ground, the event function will be zero and Diffsol will stop the integration, and we can update the velocity of the ball accordingly.
In code, the bouncing ball problem can be solved using Diffsol as follows:
use diffsol::{CraneliftJitModule, OdeBuilder, OdeSolverMethod, OdeSolverStopReason, Vector};
use plotly::{common::Mode, layout::Axis, layout::Layout, Plot, Scatter};
use std::fs;
type M = diffsol::NalgebraMat<f64>;
type CG = CraneliftJitModule;
type LS = diffsol::NalgebraLU<f64>;
fn main() {
let e = 0.8;
let problem = OdeBuilder::<M>::new()
.build_from_diffsl::<CG>(
"
g { 9.81 } h { 10.0 }
u_i {
position = h,
velocity = 0,
}
F_i {
velocity,
-g,
}
stop {
position,
}
",
)
.unwrap();
let mut solver = problem.bdf::<LS>().unwrap();
let mut position = Vec::new();
let mut velocity = Vec::new();
let mut t = Vec::new();
let final_time = 10.0;
// save the initial state
position.push(solver.state().y[0]);
velocity.push(solver.state().y[1]);
t.push(0.0);
// solve until the final time is reached
solver.set_stop_time(final_time).unwrap();
loop {
match solver.step() {
Ok(OdeSolverStopReason::InternalTimestep) => (),
Ok(OdeSolverStopReason::RootFound(t, _idx)) => {
// get the state when the event occurred
let mut y = solver.interpolate(t).unwrap();
// update the velocity of the ball
y[1] *= -e;
// make sure the ball is above the ground
y[0] = y[0].max(f64::EPSILON);
// set the state to the updated state
solver.state_mut().y.copy_from(&y);
solver.state_mut().dy[0] = y[1];
*solver.state_mut().t = t;
}
Ok(OdeSolverStopReason::TstopReached) => break,
Err(_) => panic!("unexpected solver error"),
}
position.push(solver.state().y[0]);
velocity.push(solver.state().y[1]);
t.push(solver.state().t);
}
let mut plot = Plot::new();
let position = Scatter::new(t.clone(), position)
.mode(Mode::Lines)
.name("position");
let velocity = Scatter::new(t, velocity).mode(Mode::Lines).name("velocity");
plot.add_trace(position);
plot.add_trace(velocity);
let layout = Layout::new()
.x_axis(Axis::new().title("t"))
.y_axis(Axis::new());
plot.set_layout(layout);
let plot_html = plot.to_inline_html(Some("bouncing-ball"));
fs::write("../src/primer/images/bouncing-ball.html", plot_html).expect("Unable to write file");
}Declarative approach
In DiffSL, we can move the event definition into the model itself. The stop tensor declares that an event occurs when x = 0, and the reset tensor declares the post-impact state. In this version we also reset the position to a very small positive value so that the next step starts just above the ground and does not immediately retrigger the same event.
Because the model provides both stop and reset, the high-level solve call now handles each bounce automatically. The example can therefore solve straight to the final time and read the returned trajectory directly.
use diffsol::{CraneliftJitModule, MatrixCommon, OdeBuilder, OdeSolverMethod};
use plotly::{common::Mode, layout::Axis, layout::Layout, Plot, Scatter};
use std::{fs, path::PathBuf};
type M = diffsol::NalgebraMat<f64>;
type CG = CraneliftJitModule;
type LS = diffsol::NalgebraLU<f64>;
fn main() {
let model = "
restitution { 0.8 } xeps { 1e-12 }
g { 9.81 } h { 10.0 }
u_i {
position = h,
velocity = 0,
}
F_i {
velocity,
-g,
}
stop_i {
position,
}
reset_i {
xeps,
-restitution * velocity,
}
";
let problem = OdeBuilder::<M>::new()
.build_from_diffsl::<CG>(model)
.unwrap();
let mut solver = problem.bdf::<LS>().unwrap();
let final_time = 10.0;
let (ys, t, _stop_reason) = solver.solve(final_time).unwrap();
let x: Vec<_> = ys.inner().row(0).into_iter().copied().collect();
let v: Vec<_> = ys.inner().row(1).into_iter().copied().collect();
let mut plot = Plot::new();
plot.add_trace(Scatter::new(t.clone(), x).mode(Mode::Lines).name("x"));
plot.add_trace(Scatter::new(t, v).mode(Mode::Lines).name("v"));
plot.set_layout(
Layout::new()
.x_axis(Axis::new().title("t"))
.y_axis(Axis::new()),
);
let output_path = PathBuf::from(env!("CARGO_MANIFEST_DIR"))
.join("../../book/src/primer/images/bouncing-ball-declarative.html");
let plot_html = plot.to_inline_html(Some("bouncing-ball-declarative"));
fs::write(output_path, plot_html).expect("Unable to write file");
}