Delay Differential Equations

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

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:

For this problem, we note that 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:

using DifferentialEquations
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

bc_model (generic function with 1 method)

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)

h (generic function with 1 method)

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

tau = 1
lags = [tau]

1-element Vector{Int64}:
 1

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)

DDEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 3-element Vector{Float64}:
 1.0
 1.0
 1.0

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())

MethodOfSteps{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}(Tsit5(; stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false),), NLFunctional{Rational{Int64}, Rational{Int64}}(1//100, 1//5, 10))

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 a 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)

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 40-element Vector{Float64}:
  0.0
  0.05836370738671701
  0.12916108569539597
  0.21568702493539388
  0.31476012203954973
  0.42647876001237145
  0.5486009950019299
  0.6801498148290146
  0.8195261166349427
  0.9657507075079192
  ⋮
  5.637340314503988
  6.0
  6.465321265034228
  6.9807287937908855
  7.578753298968584
  8.241987232937237
  8.98383276855844
  9.794183131169206
 10.0
u: 40-element Vector{Vector{Float64}}:
 [1.0, 1.0, 1.0]
 [0.7447276972023367, 0.9674847415959038, 0.9739922644849844]
 [0.520865690149123, 0.9216176741982295, 0.9429345997713768]
 [0.33647961749786465, 0.860710765286255, 0.9053954304799672]
 [0.20402080446565815, 0.7893156306081849, 0.8627636946755286]
 [0.11605349292143427, 0.710940127223295, 0.8151154583589272]
 [0.06263661020799872, 0.6308046470869514, 0.7637408792490028]
 [0.03223510074296827, 0.5524928215361826, 0.7096033440842167]
 [0.01594663432245627, 0.47891242866857714, 0.654067407727854]
 [0.007620742290428411, 0.41158157168839415, 0.5982822197402151]
 ⋮
 [1.9479598719208864e-12, 0.0025766300373552546, 0.0197374598472333]
 [3.3296189152665293e-13, 0.001729122637428431, 0.014435113877197158]
 [5.5064416255269584e-14, 0.0010364504244003874, 0.009606890211128613]
 [1.1838460929002524e-14, 0.0005879461631044232, 0.006079732552649997]
 [5.109127121132652e-15, 0.0003045520635702376, 0.003549648850830764]
 [4.075064652637717e-15, 0.0001468398691910801, 0.0019394744851741882]
 [6.623872359855281e-15, 6.493848650511485e-5, 0.0009786579729539678]
 [1.912280041308691e-14, 2.6637039533047875e-5, 0.0004599558275824244]
 [6.69983491663282e-15, 2.1240382256782126e-5, 0.0003792544387315653]

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)

h (generic function with 2 methods)

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

bc_model (generic function with 1 method)

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

bc_model (generic function with 1 method)

Note that this requires that we define the historical values

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

h (generic function with 2 methods)

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)

h (generic function with 3 methods)

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 also discussed 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.

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)

retcode: Success
Interpolation: 3rd order Hermite
t: 56-element Vector{Float64}:
  0.0
  0.028686567733943524
  0.05496370797935472
  0.09262168203274358
  0.1328816015643946
  0.1804089117440632
  0.23156096928227454
  0.28745950203007636
  0.34649352838806813
  0.40865997317696995
  ⋮
  6.709304549678902
  7.077771808058965
  7.457378129664066
  7.849229500378771
  8.253885262007339
  8.671507191596193
  9.102489930712723
  9.548109990519311
 10.0
u: 56-element Vector{Vector{Float64}}:
 [1.0, 1.0, 1.0]
 [0.8651377233256287, 0.9847919107281841, 0.9871577813957303]
 [0.7576257549925667, 0.9695375729958094, 0.9754952861529256]
 [0.6264188253434527, 0.9459514608421642, 0.9589125735972589]
 [0.511173996718998, 0.9190777328442522, 0.9413127662016809]
 [0.4021017755857873, 0.8859139628064613, 0.9206598661690842]
 [0.31056686797828226, 0.8492819001921053, 0.8985424920191215]
 [0.23418948608094928, 0.8089369051538083, 0.8744786633816477]
 [0.17382236795861314, 0.766681887677747, 0.8491794697919689]
 [0.1269922932046665, 0.7231641434076941, 0.8226801184713568]
 ⋮
 [7.633678240249587e-14, 0.0007930624927795012, 0.007743143315919018]
 [2.290654525439246e-14, 0.00052886086662893, 0.005575060518185558]
 [7.192689065388156e-15, 0.000348386861669049, 0.00396242735079288]
 [2.3938297589461707e-15, 0.00022643491366957353, 0.0027774610777958984]
 [8.537747837490129e-16, 0.00014511825271079198, 0.001919152324510346]
 [3.290737503207264e-16, 9.169024318664183e-5, 0.0013070717660449717]
 [1.3830954331558912e-16, 5.7090180277528864e-5, 0.0008771566981009424]
 [6.431419123095559e-17, 3.498102566610629e-5, 0.0005793516695673948]
 [3.127939839039937e-17, 2.128729410362161e-5, 0.00037958433259645345]

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 passing 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)

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 31-element Vector{Float64}:
  0.0
  0.05836370738671701
  0.12916108569539597
  0.21568702493539388
  0.31476012203954973
  0.42647876001237145
  0.5486009950019299
  0.6801498148290146
  0.8195261166349427
  0.9657507075079192
  ⋮
  4.70131900781237
  5.319727380784704
  5.975944439730393
  6.672902478790312
  7.412222703550066
  8.203867606081308
  9.055215255977002
  9.965347765311543
 10.0
u: 31-element Vector{Vector{Float64}}:
 [1.0, 1.0, 1.0]
 [0.7447276972023367, 0.9674847415959038, 0.9739922644849844]
 [0.520865690149123, 0.9216176741982295, 0.9429345997713768]
 [0.33647961749786465, 0.860710765286255, 0.9053954304799672]
 [0.20402080446565815, 0.7893156306081849, 0.8627636946755286]
 [0.11605349292143427, 0.710940127223295, 0.8151154583589272]
 [0.06263661020799872, 0.6308046470869514, 0.7637408792490028]
 [0.03223510074296827, 0.5524928215361826, 0.7096033440842167]
 [0.01594663432245627, 0.47891242866857714, 0.654067407727854]
 [0.007620742290428411, 0.41158157168839415, 0.5982822197402151]
 ⋮
 [1.637789748146323e-9, 0.007211518284286898, 0.043270644077427596]
 [8.555301578118273e-10, 0.0036540032665950783, 0.025862653910404754]
 [6.387919765560543e-10, 0.0017756654863386229, 0.014739430025735469]
 [6.957421403133639e-10, 0.000825017766449486, 0.007996169982628906]
 [1.1062342089489479e-9, 0.0003658747495832752, 0.004126394925649491]
 [2.7421619004585237e-9, 0.0001531915519910343, 0.0020084764604955713]
 [1.0932031990743894e-8, 6.0070746916231975e-5, 0.0009160350421131426]
 [6.748471780530628e-8, 2.2069323013970547e-5, 0.00039185666579902566]
 [5.655280359438763e-8, 2.1246239705310403e-5, 0.00037932174234340285]

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.