Simulate a bouncing ball?

Question

Is it possible to create a simple model of a bouncing ball, using Julia's equation solvers?

I started with this:

using ODE

function bb(t, f)
    (y, v) = f
    dy_dt = v
    dv_dt = -9.81
    [dy_dt, dv_dt]
end

const y0 =  50.0             # height
const v0 =   0.0             # velocity
const startpos = [y0; v0]

ts = 0.0:0.25:10             # time span

t, res = ode45(bb, startpos, ts)

which produces useful-looking numbers:

julia> t
44-element Array{Float64,1}:
  0.0
  0.0551392
  0.25
  0.5
  0.75
  1.0
  ⋮
  8.75
  9.0
  9.25
  9.5
  9.75
 10.0

julia> res
44-element Array{Array{Float64,1},1}:
 [50.0,0.0]
 [49.9851,-0.540915]
 [49.6934,-2.4525]
 [48.7738,-4.905]
 [47.2409,-7.3575]
 ⋮
 [-392.676,-93.195]
 [-416.282,-95.6475]
 [-440.5,-98.1]

But somehow it needs to intervene when the height is 0, and reverse the velocity. Or am I on the wrong track?

Answer 1

DifferentialEquations.jl offers sophisticated callbacks and event handling. Since the DifferentialEquations.jl algorithms are about 10x faster while offering a higher order interpolation, these algorithms are clearly the better choose here anyways.

The first link is the documentation which shows how to do the event handling. The easy interface uses the macros. I start by defining the function.

f = @ode_def BallBounce begin
  dy =  v
  dv = -g
end g=9.81

Here I am showing ParameterizedFunctions.jl to make the syntax nicer, but you can define the function directly as an in-place update f(t,u,du) (like Sundials.jl). Next you define the function which determines when an event takes place. It can be any function which is positive and hits zero at the event time. Here, we are checking for when the ball hits the ground, or for when y=0, so:

function event_f(t,u) # Event when event_f(t,u,k) == 0
  u[1]
end

Next you say what to do when the event occurs. Here we want to reverse the sign of the velocity:

function apply_event!(u,cache)
  u[2] = -u[2]
end

You put these functions together to build the callback using the macros:

callback = @ode_callback begin
  @ode_event event_f apply_event!
end

Now you solve as usual. You define the ODEProblem using f and the initial condition, and you call solve on a timespan. The only thing extra is you pass the callback along with the solver:

u0 = [50.0,0.0]
prob = ODEProblem(f,u0)
tspan = [0;15]
sol = solve(prob,tspan,callback=callback)

Then we can use the plot recipe to automatically plot the solution:

plot(sol)

The result is this:

A few things to notice here:

1. DifferentialEquations.jl will automatically use an interpolation to more safely check for the event. For example, if the event happened within a timestep but not at the ends, DifferentialEquations.jl will still find it. More or less interpolations points can be included as options to the @ode_event macro.

2. DifferentialEquations.jl used a rootfinding method to hone in on the moment of the event. Even though the adaptive solver steps past the event, by using rootfinding on the interpolation it finds the exact time of the event, and thus gets the discontinuity right. You can see that in the graph since the ball never goes negative.

3. There is a whole lot more this can do. Check out the docs. You can do pretty much anything with this. For example, have your ODE changing size over the run to model a population of cells with birth and deaths. This is something other solver packages can't do.

4. Even with all of these features, speed is not compromised.

Let me know if you need any extra functionality added to the "ease of use" interface macros.