# Bouncing Ball

Another basic dynamics model is the bouncing ball model. This introduces a few interesting features of the language. First, this will involve a few more variables and the model itself is a second order system. But in addition, such a model must handle discrete events (the ball bouncing) and reinitialization of state variables (specifically), the ball's vertical velocity.

Let's start with the basic physics which are captured by the following system of equations:

using ModelingToolkit
using DifferentialEquations
using Plots

@variables t h(t) v(t)
D = Differential(t)

@parameters M=1.0

eqs = [
    M*D(v) ~ -M*9.81
    v ~ D(h)
]

@named sys = ODESystem(eqs, t)

prob = ODEProblem(structural_simplify(sys), [h => 10.0, v => 0], (0, 10), [], jac = true)

sol = solve(prob)
display(plot(sol, idxs=[h]))

One difference from our previous example is that now we have multiple equations. Even in our previous example we used the Julia vector notation but it is easy to miss that it was, in fact, a vector since there was only one equation in it. In this case, however, we see that there are two equations where each equation is simply a component in the eqs vector.

Note that with this model, there is no surface to bounce off of. So the ball just continues to fall, i.e.,

Now, let's consider the following augmented model with a few minor changes (highlighted):

using ModelingToolkit
using DifferentialEquations
using Plots

@variables t h(t) v(t)
D = Differential(t)

@parameters M=1.0

when = [h ~ 0] 
what = [v ~ -0.9*v]

eqs = [
    M*D(v) ~ -M*9.81
    v ~ D(h)
]

@named sys = ODESystem(eqs, t; continuous_events = when => what)

prob = ODEProblem(structural_simplify(sys), [h => 10.0, v => 0], (0, 10), [], jac = true)

sol = solve(prob)
display(plot(sol, idxs=[h]))

The when variable holds a vector of conditions. When those conditions become true, an event will be generated. What the event does, when it occurs, is indicated by the what variable. This is indicating how we can reinitialize v. So in this case, v will take on the value of -0.9 times the previous value. The 0.9 in this case is the so-called "coefficient of restitution" and indicates how much of the ball's momentum is retained by the ball after it bounces.

With this change, the ball actually starts bouncing:

Note that if we run this model out, we get "Zeno" behavior (events happening at higher and higher frequencies) which makes it very difficult for the solver to make progress. So for now, don't try to simulate out much further in time.