📂Julia

How to Terminate Differential Equation Solving Early in Julia

Overview

SIR model: $$ \begin{align*} {{d S} \over {d t}} =& - {{ \beta } \over { N }} I S \\ {{d I} \over {d t}} =& {{ \beta } \over { N }} S I - \mu I \\ {{d R} \over {d t}} =& \mu I \end{align*} $$

For instance, when simulating an infectious disease spread model like the one above, the fixed point will be at $I = 0$, and if we consider the number of infected people to have gotten sufficiently close to $0$, we might want to terminate the simulation at that point.

Code

In Julia’s differential equation solving package, among the features of the event handler¹, one can use terminate! to terminate the solving of a differential equation early².

Existing

prob = ODEProblem(SIR, rand(3), (0, 200))
sol = solve(prob)

If we solve the differential equation over the interval $t \in [0, 200]$ from an arbitrary initial condition, the solution may look like the above. As can be seen, there is no particularly meaningful change after $t > 50$, and if such a simulation is repeated countless times, there is no need to solve all the way to $t = 200$.

Early Termination

condition(u,t,integrator) = u[2] - 0.01
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition,affect!)
sol = solve(prob, callback = cb);

On the other hand, when there is a callback instructing to terminate under the condition $I < 0.01$, affect! is called at the point where the sign of the condition function changes, and the differential equation solving is terminated by the terminate! called within it.

Complete Code

using DifferentialEquations, Plots

function SIR(du, u, p, t)
    S, I, R = u
    du[1] = -0.5 * S * I
    du[2] = 0.5 * S * I - 0.1 * I
    du[3] = 0.1 * I
end
prob = ODEProblem(SIR, rand(3), (0, 200))

sol = solve(prob)
plot(sol)

condition(u,t,integrator) = u[2] - 0.01
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition,affect!)
sol = solve(prob, callback = cb);
plot(sol)