Delay Differential Equations

This tutorial will introduce you to the functionality for solving delay differential equations.

Note

This tutorial assumes you have read the Ordinary Differential Equations tutorial.

Delay differential equations are equations which have a delayed argument. To allow for specifying the delayed argument, the function definition for a delay differential equation is expanded to include a history function h(p, t) which uses interpolations throughout the solution's history to form a continuous extension of the solver's past and depends on parameters p and time t. The function signature for a delay differential equation is f(u, h, p, t) for not in-place computations, and f(du, u, h, p, t) for in-place computations.

In this example we will solve a model of breast cancer growth kinetics:

\begin{aligned} dx_{0} &= \frac{v_{0}}{1+\beta_{0}\left(x_{2}(t-\tau)\right)^{2}}\left(p_{0}-q_{0}\right)x_{0}(t)-d_{0}x_{0}(t)\\ dx_{1} &= \frac{v_{0}}{1+\beta_{0}\left(x_{2}(t-\tau)\right)^{2}}\left(1-p_{0}+q_{0}\right)x_{0}(t)\\ &+ \frac{v_{1}}{1+\beta_{1}\left(x_{2}(t-\tau)\right)^{2}}\left(p_{1}-q_{1}\right)x_{1}(t)-d_{1}x_{1}(t)\\ dx_{2} &= \frac{v_{1}}{1+\beta_{1}\left(x_{2}(t-\tau)\right)^{2}}\left(1-p_{1}+q_{1}\right)x_{1}(t)-d_{2}x_{2}(t) \end{aligned}

For this problem we note that $\tau$ is constant, and thus we can use a method which exploits this behavior. We first write out the equation using the appropriate function signature. Most of the equation writing is the same, though we use the history function by first interpolating and then choosing the components. Thus the ith component at time t-tau is given by h(p, t-tau)[i]. Components with no delays are written as in the ODE.

Thus, the function for this model is given by:

function bc_model(du,u,h,p,t)
p0,q0,v0,d0,p1,q1,v1,d1,d2,beta0,beta1,tau = p
hist3 = h(p, t-tau)[3]
du[1] = (v0/(1+beta0*(hist3^2))) * (p0 - q0)*u[1] - d0*u[1]
du[2] = (v0/(1+beta0*(hist3^2))) * (1 - p0 + q0)*u[1] +
(v1/(1+beta1*(hist3^2))) * (p1 - q1)*u[2] - d1*u[2]
du[3] = (v1/(1+beta1*(hist3^2))) * (1 - p1 + q1)*u[2] - d2*u[3]
end

Now we build a DDEProblem. The signature

prob = DDEProblem(f, u0, h, tspan, p=SciMLBase.NullParameters();
constant_lags=[], dependent_lags=[], kwargs...)

is very similar to ODEs, where we now have to give the lags and a function h. h is the history function that declares what the values were before the time the model starts. Here we will assume that for all time before t0 the values were 1 and define h as an out-of-place function:

h(p, t) = ones(3)

To use the constant lag model, we have to declare the lags. Here we will use tau=1.

tau = 1
lags = [tau]

Next, we choose to solve on the timespan (0.0,10.0) and create the problem type:

p0 = 0.2; q0 = 0.3; v0 = 1; d0 = 5
p1 = 0.2; q1 = 0.3; v1 = 1; d1 = 1
d2 = 1; beta0 = 1; beta1 = 1
p = (p0,q0,v0,d0,p1,q1,v1,d1,d2,beta0,beta1,tau)
tspan = (0.0,10.0)
u0 = [1.0,1.0,1.0]

prob = DDEProblem(bc_model,u0,h,tspan,p; constant_lags=lags)

An efficient way to solve this problem (given the constant lags) is with the MethodOfSteps solver. Through the magic that is Julia, it translates an OrdinaryDiffEq.jl ODE solver method into a method for delay differential equations which is highly efficient due to sweet compiler magic. A good choice is the order 5 method Tsit5():

alg = MethodOfSteps(Tsit5())

For lower tolerance solving, one can use the BS3() algorithm to good effect (this combination is similar to the MATLAB dde23, but more efficient tableau), and for high tolerances the Vern6() algorithm will give an 6th order solution.

To solve the problem with this algorithm, we do the same thing we'd do with other methods on the common interface:

sol = solve(prob,alg)

Note that everything available to OrdinaryDiffEq.jl can be used here, including event handling and other callbacks. The solution object has the same interface as for ODEs. For example, we can use the same plot recipes to view the results:

using Plots; plot(sol)

Speeding Up Interpolations with Idxs

We can speed up the previous problem in two different ways. First of all, if we need to interpolate multiple values from a previous time, we can use the in-place form for the history function h(out, p, t) which writes the output to out. In this case, we must supply the history initial conditions as in-place as well. For the previous example, that's simply

h(out, p, t) = (out.=1.0)

and then our DDE is:

const out = zeros(3) # Define a cache variable
function bc_model(du,u,h,p,t)
h(out, p, t-tau) # updates out to be the correct history function
du[1] = (v0/(1+beta0*(out[3]^2))) * (p0 - q0)*u[1] - d0*u[1]
du[2] = (v0/(1+beta0*(out[3]^2))) * (1 - p0 + q0)*u[1] +
(v1/(1+beta1*(out[3]^2))) * (p1 - q1)*u[2] - d1*u[2]
du[3] = (v1/(1+beta1*(out[3]^2))) * (1 - p1 + q1)*u[2] - d2*u[3]
end

However, we can do something even slicker in most cases. We only ever needed to interpolate past values at index 3. Instead of generating a bunch of arrays, we can instead ask specifically for that value by passing the keyword idxs = 3. The DDE function bc_model is now:

function bc_model(du,u,h,p,t)
u3_past_sq = h(p, t-tau; idxs=3)^2
du[1] = (v0/(1+beta0*(u3_past_sq))) * (p0 - q0)*u[1] - d0*u[1]
du[2] = (v0/(1+beta0*(u3_past_sq))) * (1 - p0 + q0)*u[1] +
(v1/(1+beta1*(u3_past_sq))) * (p1 - q1)*u[2] - d1*u[2]
du[3] = (v1/(1+beta1*(u3_past_sq))) * (1 - p1 + q1)*u[2] - d2*u[3]
end

Note that this requires that we define the historical values

h(p, t; idxs=nothing) = typeof(idxs) <: Number ? 1.0 : ones(3)

where idxs can be an integer for which variable in the history to compute, and here for any number idxs we give back 1.0. Note that if we wanted to use past values of the ith derivative then we would call the history function h(p, t, Val{i}) in our DDE function and would have to define a dispatch like

h(p, t, ::Type{Val{1}}) = zeros(3)

to say that derivatives before t0 are zero for any index. Again, we could use an in-place function instead or only compute specific indices by passing an idxs keyword.

The functional forms for the history function are discussed also on the DDEProblem page.

Undeclared Delays and State-Dependent Delays via Residual Control

You might have noticed DifferentialEquations.jl allows you to solve problems with undeclared delays since you can interpolate h at any value. This is a feature, but use it with caution. Undeclared delays can increase the error in the solution. It's recommended that you use a method with a residual control, such as MethodOfSteps(RK4()) whenever there are undeclared delays. With this you can use interpolated derivatives, solve functional differential equations by using quadrature on the interpolant, etc. However, note that residual control solves with a low level of accuracy, so the tolerances should be made very small and the solution should not be trusted for more than 2-3 decimal places.

Note: MethodOfSteps(RK4()) with undeclared delays is similar to MATLAB's ddesd. Thus, for example, the following is similar to solving the example from above with residual control:

prob = DDEProblem(bc_model,u0,h,tspan)
alg = MethodOfSteps(RK4())
sol = solve(prob,alg)

Note that this method can solve problems with state-dependent delays.

State-Dependent Delay Discontinuity Tracking

State-dependent delays are problems where the delay is allowed to be a function of the current state. They can be more efficiently solved with discontinuity tracking. To do this in DifferentialEquations.jl, requires to pass lag functions g(u,p,t) as keyword dependent_lags to the DDEProblem definition. Other than that, everything else is the same, and one solves that problem using the common interface.

We can solve the above problem with dependent delay tracking by declaring the dependent lags and solving with a MethodOfSteps algorithm:

prob = DDEProblem(bc_model,u0,h,tspan; dependent_lags = ((u,p,t) -> tau,))
alg = MethodOfSteps(Tsit5())
sol = solve(prob,alg)

Here we treated the single lag t-tau as a state-dependent delay. Of course, you can then replace that tuple of functions with whatever functions match your lags.