In working with a differential equation, our system will evolve through many states. Particular states of the system may be of interest to us, and we say that an ***"event"*** is triggered when our system reaches these states. For example, events may include the moment when our system reaches a particular temperature or velocity. We ***handle*** these events with ***callbacks***, which tell us what to do once an event has been triggered.

These callbacks allow for a lot more than event handling, however. For example, we can use callbacks to achieve high-level behavior like exactly preserve conservation laws and save the trace of a matrix at pre-defined time points. This extra functionality allows us to use the callback system as a modding system for the DiffEq ecosystem's solvers.

This tutorial is an introduction to the callback and event handling system in DifferentialEquations.jl, documented in the Event Handling and Callback Functions page of the documentation. We will also introduce you to some of the most widely used callbacks in the Callback Library, which is a library of pre-built mods.

Event handling is done through continuous callbacks. Callbacks take a function, `condition`

, which triggers an `affect!`

when `condition == 0`

. These callbacks are called "continuous" because they will utilize rootfinding on the interpolation to find the "exact" time point at which the condition takes place and apply the `affect!`

at that time point.

***Let's use a bouncing ball as a simple system to explain events and callbacks.*** Let's take Newton's model of a ball falling towards the Earth's surface via a gravitational constant `g`

. In this case, the velocity is changing via `-g`

, and position is changing via the velocity. Therefore we receive the system of ODEs:

using DifferentialEquations, ParameterizedFunctions ball! = @ode_def BallBounce begin dy = v dv = -g end g

(::Main.WeaveSandBox8.BallBounce{getfield(Main.WeaveSandBox8, Symbol("##1#5 ")),getfield(Main.WeaveSandBox8, Symbol("##2#6")),getfield(Main.WeaveSandBo x8, Symbol("##3#7")),Nothing,Nothing,getfield(Main.WeaveSandBox8, Symbol("# #4#8")),Expr,Expr}) (generic function with 2 methods)

We want the callback to trigger when `y=0`

since that's when the ball will hit the Earth's surface (our event). We do this with the condition:

function condition(u,t,integrator) u[1] end

condition (generic function with 1 method)

Recall that the `condition`

will trigger when it evaluates to zero, and here it will evaluate to zero when `u[1] == 0`

, which occurs when `v == 0`

. *Now we have to say what we want the callback to do.* Callbacks make use of the Integrator Interface. Instead of giving a full description, a quick and usable rundown is:

Values are strored in

`integrator.u`

Times are stored in

`integrator.t`

The parameters are stored in

`integrator.p`

`integrator(t)`

performs an interpolation in the current interval between`integrator.tprev`

and`integrator.t`

(and allows extrapolation)User-defined options (tolerances, etc.) are stored in

`integrator.opts`

`integrator.sol`

is the current solution object. Note that`integrator.sol.prob`

is the current problem

While there's a lot more on the integrator interface page, that's a working knowledge of what's there.

What we want to do with our `affect!`

is to "make the ball bounce". Mathematically speaking, the ball bounces when the sign of the velocity flips. As an added behavior, let's also use a small friction constant to dampen the ball's velocity. This way only a percentage of the velocity will be retained when the event is triggered and the callback is used. We'll define this behavior in the `affect!`

function:

function affect!(integrator) integrator.u[2] = -integrator.p[2] * integrator.u[2] end

affect! (generic function with 1 method)

`integrator.u[2]`

is the second value of our model, which is `v`

or velocity, and `integrator.p[2]`

, is our friction coefficient.

Therefore `affect!`

can be read as follows: `affect!`

will take the current value of velocity, and multiply it `-1`

multiplied by our friction coefficient. Therefore the ball will change direction and its velocity will dampen when `affect!`

is called.

Now let's build the `ContinuousCallback`

:

bounce_cb = ContinuousCallback(condition,affect!)

DiffEqBase.ContinuousCallback{typeof(Main.WeaveSandBox8.condition),typeof(M ain.WeaveSandBox8.affect!),typeof(Main.WeaveSandBox8.affect!),typeof(DiffEq Base.INITIALIZE_DEFAULT),Float64,Int64,Nothing}(Main.WeaveSandBox8.conditio n, Main.WeaveSandBox8.affect!, Main.WeaveSandBox8.affect!, DiffEqBase.INITI ALIZE_DEFAULT, nothing, true, 10, Bool[true, true], 2.220446049250313e-15, 0)

Now let's make an `ODEProblem`

which has our callback:

u0 = [50.0,0.0] tspan = (0.0,15.0) p = (9.8,0.9) prob = ODEProblem(ball!,u0,tspan,p,callback=bounce_cb)

ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true timespan: (0.0, 15.0) u0: [50.0, 0.0]

Notice that we chose a friction constant of `0.9`

. Now we can solve the problem and plot the solution as we normally would:

sol = solve(prob,Tsit5()) using Plots; gr() plot(sol)

and tada, the ball bounces! Notice that the `ContinuousCallback`

is using the interpolation to apply the effect "exactly" when `v == 0`

. This is crucial for model correctness, and thus when this property is needed a `ContinuousCallback`

should be used.

In our example we used a constant coefficient of friction, but if we are bouncing the ball in the same place we may be smoothing the surface (say, squishing the grass), causing there to be less friction after each bounce. In this more advanced model, we want the friction coefficient at the next bounce to be `sqrt(friction)`

from the previous bounce (since `friction < 1`

, `sqrt(friction) > friction`

and `sqrt(friction) < 1`

).

Hint: there are many ways to implement this. One way to do it is to make `p`

a `Vector`

and mutate the friction coefficient in the `affect!`

.

A discrete callback checks a `condition`

after every integration step and, if true, it will apply an `affect!`

. For example, let's say that at time `t=2`

we want to include that a kid kicked the ball, adding `20`

to the current velocity. This kind of situation, where we want to add a specific behavior which does not require rootfinding, is a good candidate for a `DiscreteCallback`

. In this case, the `condition`

is a boolean for whether to apply the `affect!`

, so:

function condition_kick(u,t,integrator) t == 2 end

condition_kick (generic function with 1 method)

We want the kick to occur at `t=2`

, so we check for that time point. When we are at this time point, we want to do:

function affect_kick!(integrator) integrator.u[2] += 50 end

affect_kick! (generic function with 1 method)

Now we build the problem as before:

kick_cb = DiscreteCallback(condition_kick,affect_kick!) u0 = [50.0,0.0] tspan = (0.0,10.0) p = (9.8,0.9) prob = ODEProblem(ball!,u0,tspan,p,callback=kick_cb)

ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true timespan: (0.0, 10.0) u0: [50.0, 0.0]

Note that, since we are requiring our effect at exactly the time `t=2`

, we need to tell the integration scheme to step at exactly `t=2`

to apply this callback. This is done via the option `tstops`

, which is like `saveat`

but means "stop at these values".

sol = solve(prob,Tsit5(),tstops=[2.0]) plot(sol)

Note that this example could've been done with a `ContinuousCallback`

by checking the condition `t-2`

.

In some cases you may want to merge callbacks to build up more complex behavior. In our previous result, notice that the model is unphysical because the ball goes below zero! What we really need to do is add the bounce callback together with the kick. This can be achieved through the `CallbackSet`

.

cb = CallbackSet(bounce_cb,kick_cb)

DiffEqBase.CallbackSet{Tuple{DiffEqBase.ContinuousCallback{typeof(Main.Weav eSandBox8.condition),typeof(Main.WeaveSandBox8.affect!),typeof(Main.WeaveSa ndBox8.affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing }},Tuple{DiffEqBase.DiscreteCallback{typeof(Main.WeaveSandBox8.condition_ki ck),typeof(Main.WeaveSandBox8.affect_kick!),typeof(DiffEqBase.INITIALIZE_DE FAULT)}}}((DiffEqBase.ContinuousCallback{typeof(Main.WeaveSandBox8.conditio n),typeof(Main.WeaveSandBox8.affect!),typeof(Main.WeaveSandBox8.affect!),ty peof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}(Main.WeaveSandBo x8.condition, Main.WeaveSandBox8.affect!, Main.WeaveSandBox8.affect!, DiffE qBase.INITIALIZE_DEFAULT, nothing, true, 10, Bool[true, true], 2.2204460492 50313e-15, 0),), (DiffEqBase.DiscreteCallback{typeof(Main.WeaveSandBox8.con dition_kick),typeof(Main.WeaveSandBox8.affect_kick!),typeof(DiffEqBase.INIT IALIZE_DEFAULT)}(Main.WeaveSandBox8.condition_kick, Main.WeaveSandBox8.affe ct_kick!, DiffEqBase.INITIALIZE_DEFAULT, Bool[true, true]),))

A `CallbackSet`

merges their behavior together. The logic is as follows. In a given interval, if there are multiple continuous callbacks that would trigger, only the one that triggers at the earliest time is used. The time is pulled back to where that continuous callback is triggered, and then the `DiscreteCallback`

s in the callback set are called in order.

u0 = [50.0,0.0] tspan = (0.0,15.0) p = (9.8,0.9) prob = ODEProblem(ball!,u0,tspan,p,callback=cb) sol = solve(prob,Tsit5(),tstops=[2.0]) plot(sol)

Notice that we have now merged the behaviors. We can then nest this as deep as we like.

Add to the model a linear wind with resistance that changes the acceleration to `-g + k*v`

after `t=10`

. Do so by adding another parameter and allowing it to be zero until a specific time point where a third callback triggers the change.

Let's look at another model now: the model of the Harmonic Oscillator. We can write this as:

u0 = [1.,0.] harmonic! = @ode_def HarmonicOscillator begin dv = -x dx = v end tspan = (0.0,10.0) prob = ODEProblem(harmonic!,u0,tspan) sol = solve(prob) plot(sol)

Let's instead stop the integration when a condition is met. From the Integrator Interface stepping controls we see that `terminate!(integrator)`

will cause the integration to end. So our new `affect!`

is simply:

function terminate_affect!(integrator) terminate!(integrator) end

terminate_affect! (generic function with 1 method)

Let's first stop the integration when the particle moves back to `x=0`

. This means we want to use the condition:

function terminate_condition(u,t,integrator) u[2] end terminate_cb = ContinuousCallback(terminate_condition,terminate_affect!)

DiffEqBase.ContinuousCallback{typeof(Main.WeaveSandBox8.terminate_condition ),typeof(Main.WeaveSandBox8.terminate_affect!),typeof(Main.WeaveSandBox8.te rminate_affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothin g}(Main.WeaveSandBox8.terminate_condition, Main.WeaveSandBox8.terminate_aff ect!, Main.WeaveSandBox8.terminate_affect!, DiffEqBase.INITIALIZE_DEFAULT, nothing, true, 10, Bool[true, true], 2.220446049250313e-15, 0)

Note that instead of adding callbacks to the problem, we can also add them to the `solve`

command. This will automatically form a `CallbackSet`

with any problem-related callbacks and naturally allows you to distinguish between model features and integration controls.

sol = solve(prob,callback=terminate_cb) plot(sol)

Notice that the harmonic oscilator's true solution here is `sin`

and `cosine`

, and thus we would expect this return to zero to happen at `t=π`

:

sol.t[end]

3.1415902498303465

This is one way to approximate π! Lower tolerances and arbitrary precision numbers can make this more exact, but let's not look at that. Instead, what if we wanted to halt the integration after exactly one cycle? To do so we would need to ignore the first zero-crossing. Luckily in these types of scenarios there's usually a structure to the problem that can be exploited. Here, we only want to trigger the `affect!`

when crossing from positive to negative, and not when crossing from negative to positive. In other words, we want our `affect!`

to only occur on upcrossings.

If the `ContinuousCallback`

constructor is given a single `affect!`

, it will occur on both upcrossings and downcrossings. If there are two `affect!`

s given, then the first is for upcrossings and the second is for downcrossings. An `affect!`

can be ignored by using `nothing`

. Together, the "upcrossing-only" version of the effect means that the first `affect!`

is what we defined above and the second is `nothing`

. Therefore we want:

terminate_upcrossing_cb = ContinuousCallback(terminate_condition,terminate_affect!,nothing)

DiffEqBase.ContinuousCallback{typeof(Main.WeaveSandBox8.terminate_condition ),typeof(Main.WeaveSandBox8.terminate_affect!),Nothing,typeof(DiffEqBase.IN ITIALIZE_DEFAULT),Float64,Int64,Nothing}(Main.WeaveSandBox8.terminate_condi tion, Main.WeaveSandBox8.terminate_affect!, nothing, DiffEqBase.INITIALIZE_ DEFAULT, nothing, true, 10, Bool[true, true], 2.220446049250313e-15, 0)

Which gives us:

sol = solve(prob,callback=terminate_upcrossing_cb) plot(sol)

As you can see, callbacks can be very useful and through `CallbackSets`

we can merge together various behaviors. Because of this utility, there is a library of pre-built callbacks known as the Callback Library. We will walk through a few examples where these callbacks can come in handy.

One callback is the manifold projection callback. Essentially, you can define any manifold `g(sol)=0`

which the solution must live on, and cause the integration to project to that manifold after every step. As an example, let's see what happens if we naively run the harmonic oscillator for a long time:

tspan = (0.0,10000.0) prob = ODEProblem(harmonic!,u0,tspan) sol = solve(prob) gr(fmt=:png) # Make it a PNG instead of an SVG since there's a lot of points! plot(sol,vars=(1,2))

plot(sol,vars=(0,1),denseplot=false)

Notice that what's going on is that the numerical solution is drifting from the true solution over this long time scale. This is because the integrator is not conserving energy.

plot(sol.t,[u[2]^2 + u[1]^2 for u in sol.u]) # Energy ~ x^2 + v^2

Some integration techniques like symplectic integrators are designed to mitigate this issue, but instead let's tackle the problem by enforcing conservation of energy. To do so, we define our manifold as the one where energy equals 1 (since that holds in the initial condition), that is:

function g(resid,u,p,t) resid[1] = u[2]^2 + u[1]^2 - 1 resid[2] = 0 end

g (generic function with 1 method)

Here the residual measures how far from our desired energy we are, and the number of conditions matches the size of our system (we ignored the second one by making the residual 0). Thus we define a `ManifoldProjection`

callback and add that to the solver:

cb = ManifoldProjection(g) sol = solve(prob,callback=cb) plot(sol,vars=(1,2))

plot(sol,vars=(0,1),denseplot=false)

Now we have "perfect" energy conservation, where if it's ever violated too much the solution will get projected back to `energy=1`

.

u1,u2 = sol[500] u2^2 + u1^2

1.0000425845786414

While choosing different integration schemes and using lower tolerances can achieve this effect as well, this can be a nice way to enforce physical constraints and is thus used in many disciplines like molecular dynamics. Another such domain constraining callback is the `PositiveCallback()`

which can be used to enforce positivity of the variables.

The `SavingCallback`

can be used to allow for special saving behavior. Let's take a linear ODE define on a system of 1000x1000 matrices:

prob = ODEProblem((du,u,p,t)->du.=u,rand(1000,1000),(0.0,1.0))

ODEProblem with uType Array{Float64,2} and tType Float64. In-place: true timespan: (0.0, 1.0) u0: [0.620858 0.0652844 … 0.791104 0.126102; 0.941786 0.411355 … 0.0193275 0.155585; … ; 0.0920818 0.250822 … 0.966273 0.292458; 0.270047 0.335093 … 0 .338701 0.826523]

In fields like quantum mechanics you may only want to know specific properties of the solution such as the trace or the norm of the matrix. Saving all of the 1000x1000 matrices can be a costly way to get this information! Instead, we can use the `SavingCallback`

to save the `trace`

and `norm`

at specified times. To do so, we first define our `SavedValues`

cache. Our time is in terms of `Float64`

, and we want to save tuples of `Float64`

s (one for the `trace`

and one for the `norm`

), and thus we generate the cache as:

saved_values = SavedValues(Float64, Tuple{Float64,Float64})

SavedValues{tType=Float64, savevalType=Tuple{Float64,Float64}} t: Float64[] saveval: Tuple{Float64,Float64}[]

Now we define the `SavingCallback`

by giving it a function of `(u,p,t,integrator)`

that returns the values to save, and the cache:

using LinearAlgebra cb = SavingCallback((u,t,integrator)->(tr(u),norm(u)), saved_values)

DiffEqBase.DiscreteCallback{getfield(DiffEqCallbacks, Symbol("##28#29")),Di ffEqCallbacks.SavingAffect{getfield(Main.WeaveSandBox8, Symbol("##21#22")), Float64,Tuple{Float64,Float64},DataStructures.BinaryHeap{Float64,DataStruct ures.LessThan},Array{Float64,1}},typeof(DiffEqCallbacks.saving_initialize)} (getfield(DiffEqCallbacks, Symbol("##28#29"))(), DiffEqCallbacks.SavingAffe ct{getfield(Main.WeaveSandBox8, Symbol("##21#22")),Float64,Tuple{Float64,Fl oat64},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Array{Flo at64,1}}(getfield(Main.WeaveSandBox8, Symbol("##21#22"))(), SavedValues{tTy pe=Float64, savevalType=Tuple{Float64,Float64}} t: Float64[] saveval: Tuple{Float64,Float64}[], DataStructures.BinaryHeap{Float64,DataStructures. LessThan}(DataStructures.LessThan(), Float64[]), Float64[], true, true, 0), DiffEqCallbacks.saving_initialize, Bool[false, false])

Here we take `u`

and save `(tr(u),norm(u))`

. When we solve with this callback:

sol = solve(prob, Tsit5(), callback=cb, save_everystep=false, save_start=false, save_end = false) # Turn off normal saving

retcode: Success Interpolation: 1st order linear t: 0-element Array{Float64,1} u: 0-element Array{Array{Float64,2},1}

Our values are stored in our `saved_values`

variable:

saved_values.t

5-element Array{Float64,1}: 0.0 0.10012880533703399 0.3483895172587412 0.6837345412350667 1.0

saved_values.saveval

5-element Array{Tuple{Float64,Float64},1}: (521.2188816231161, 577.5708108810427) (576.1101512173013, 638.396686933491) (738.4545731192429, 818.2930849840106) (1032.671678033474, 1144.3196696910056) (1416.8197522959554, 1570.000170864033)

By default this happened only at the solver's steps. But the `SavingCallback`

has similar controls as the integrator. For example, if we want to save at every `0.1`

seconds, we do can so using `saveat`

:

saved_values = SavedValues(Float64, Tuple{Float64,Float64}) # New cache cb = SavingCallback((u,t,integrator)->(tr(u),norm(u)), saved_values, saveat = 0.0:0.1:1.0) sol = solve(prob, Tsit5(), callback=cb, save_everystep=false, save_start=false, save_end = false) # Turn off normal saving

retcode: Success Interpolation: 1st order linear t: 0-element Array{Float64,1} u: 0-element Array{Array{Float64,2},1}

saved_values.t

11-element Array{Float64,1}: 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

saved_values.saveval

11-element Array{Tuple{Float64,Float64},1}: (521.2188816231161, 577.5708108810427) (576.0359499339363, 638.3144633285659) (636.6182023672519, 705.446606649904) (703.5718556709452, 779.6389991235587) (777.5673418091517, 861.6345569155045) (859.3446126138218, 952.2532322466953) (949.7223545366644, 1052.4022244041187) (1049.6059336885553, 1163.0847837634562) (1159.9940460767452, 1285.4075810211646) (1281.9911707976896, 1420.5945067726018) (1416.8197522959554, 1570.000170864033)

Go back to the Harmonic oscillator. Use the `SavingCallback`

to save an array for the energy over time, and do this both with and without the `ManifoldProjection`

. Plot the results to see the difference the projection makes.

This tutorial is part of the DiffEqTutorials.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqTutorials.jl

To locally run this tutorial, do the following commands:

```
using DiffEqTutorials
DiffEqTutorials.weave_file("introduction","04-callbacks_and_events.jmd")
```

Computer Information:

```
Julia Version 1.1.1
Commit 55e36cc308 (2019-05-16 04:10 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-6.0.1 (ORCJIT, ivybridge)
```

Package Information:

```
Status `~/.julia/environments/v1.1/Project.toml`
[7e558dbc-694d-5a72-987c-6f4ebed21442] ArbNumerics 0.5.4
[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.4.2
[be33ccc6-a3ff-5ff2-a52e-74243cff1e17] CUDAnative 2.2.0
[3a865a2d-5b23-5a0f-bc46-62713ec82fae] CuArrays 1.0.2
[55939f99-70c6-5e9b-8bb0-5071ed7d61fd] DecFP 0.4.8
[abce61dc-4473-55a0-ba07-351d65e31d42] Decimals 0.4.0
[ebbdde9d-f333-5424-9be2-dbf1e9acfb5e] DiffEqBayes 1.1.0
[eb300fae-53e8-50a0-950c-e21f52c2b7e0] DiffEqBiological 3.8.2
[459566f4-90b8-5000-8ac3-15dfb0a30def] DiffEqCallbacks 2.5.2
[f3b72e0c-5b89-59e1-b016-84e28bfd966d] DiffEqDevTools 2.9.0
[1130ab10-4a5a-5621-a13d-e4788d82bd4c] DiffEqParamEstim 1.6.0
[055956cb-9e8b-5191-98cc-73ae4a59e68a] DiffEqPhysics 3.1.0
[6d1b261a-3be8-11e9-3f2f-0b112a9a8436] DiffEqTutorials 0.1.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.4.0
[31c24e10-a181-5473-b8eb-7969acd0382f] Distributions 0.20.0
[497a8b3b-efae-58df-a0af-a86822472b78] DoubleFloats 0.9.1
[f6369f11-7733-5829-9624-2563aa707210] ForwardDiff 0.10.3
[c91e804a-d5a3-530f-b6f0-dfbca275c004] Gadfly 1.0.1
[7073ff75-c697-5162-941a-fcdaad2a7d2a] IJulia 1.18.1
[4138dd39-2aa7-5051-a626-17a0bb65d9c8] JLD 0.9.1
[23fbe1c1-3f47-55db-b15f-69d7ec21a316] Latexify 0.8.2
[eff96d63-e80a-5855-80a2-b1b0885c5ab7] Measurements 2.0.0
[961ee093-0014-501f-94e3-6117800e7a78] ModelingToolkit 0.2.0
[76087f3c-5699-56af-9a33-bf431cd00edd] NLopt 0.5.1
[2774e3e8-f4cf-5e23-947b-6d7e65073b56] NLsolve 4.0.0
[429524aa-4258-5aef-a3af-852621145aeb] Optim 0.18.1
[1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] OrdinaryDiffEq 5.8.1
[65888b18-ceab-5e60-b2b9-181511a3b968] ParameterizedFunctions 4.1.1
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 0.25.1
[d330b81b-6aea-500a-939a-2ce795aea3ee] PyPlot 2.8.1
[731186ca-8d62-57ce-b412-fbd966d074cd] RecursiveArrayTools 0.20.0
[90137ffa-7385-5640-81b9-e52037218182] StaticArrays 0.11.0
[f3b207a7-027a-5e70-b257-86293d7955fd] StatsPlots 0.11.0
[c3572dad-4567-51f8-b174-8c6c989267f4] Sundials 3.6.1
[1986cc42-f94f-5a68-af5c-568840ba703d] Unitful 0.15.0
[44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9] Weave 0.9.0
[b77e0a4c-d291-57a0-90e8-8db25a27a240] InteractiveUtils
[37e2e46d-f89d-539d-b4ee-838fcccc9c8e] LinearAlgebra
[44cfe95a-1eb2-52ea-b672-e2afdf69b78f] Pkg
```