The Azimuth Project
Numerical methods for delay differential equations


Sometimes it is useful to know what happens under the hood of numerical methods. The aim of this page is to provide some intuition behind the main techniques to solve delay differential equations numerically.

Preliminary definitions

A delay differential equation is a kind of differential equation where the derivative of the unknown function at any time depends on the past history of the function. More formally, let 0τ 1<τ 2,,<τ m0 \leq \tau_{1} \lt \tau_{2} , \ldots , \lt \tau_{m} be a sequence of positive delays. Then a delay differential equation (DDE) is defined as follows,

(1)y(t)=f(t,y(tτ 1),,y(tτ m))y'(t) = f(t, y(t-\tau_{1}), \ldots , y(t-\tau_{m}) )
(2)y(t)=g(t),t 0τ<tt 0y(t) = g(t), \; t_{0} - \tau \lt t \leq t_{0}

where gg is given and represents the values of yy for t 0τ<tt 0t_{0} - \tau \lt t \leq t_{0}. More about the theory of delay differential equations is discussed here.

Method of steps

We consider a particular case of delay differential equations where there is only one delayed term with a constant delay τ\tau.

(3)y(t)=f(t,y(t),y(tτ))y'(t) = f(t, y(t), y(t-\tau))
(4)y(t)=g(t),t 0τ<tt 0y(t) = g(t), \; t_{0} - \tau \lt t \leq t_{0}

We want to find a numerical solution for the problem in the interval [t 0,kτ][t_{0}, k \tau] with k >0k \in \mathbb{N}_{\gt 0}.

The basic idea is to run a traditional numerical method for ODEs (e.g.\ Runge Kutta) for each interval [t 0,t 0+τ][t_{0}, t_{0}+\tau], [t 0+τ,t 0+2τ][t_{0} + \tau, t_{0}+2 \tau] and so on. In fact, when we solve the problem for [t 0+iτ,t 0+(i+1)τ][t_{0} + i \tau, t_{0}+(i+1) \tau], we already know the values of y(t)y(t) for tt 0+iτt \leq t_{0} + i \tau from previous steps.

The following code illustrates the basic algorithm in pseudo-Python.

# initial conditions
y0 = 0.55
# the number of intervals
k = 2
# the number of points in an interval
n = 5
# the delay
tau = 1.
# set the step size
h = tau/n
# the history function for t in [t0-tau, t0]
# (for simplicity's sake we take a constant function)
hist = lambda t: y0
# the global solution for [t0-tau, k*tau] is stored in an array
# (initialization with values of the history function)
sol = array( [hist(i*h) for i in range( n*(k+1) + 1 )] )

# for each interval solve an IVP using your favorite solver
for i in range(k-1):
     # p is the function p(t) = y(t-tau)
     p = lambda t:sol[  floor( n*(1 + (t-tau)/tau) ) ]   
     # solve the following ivp 
     let lsol be the local solution for 
        y'(t) = f(t, y(t), p(t)), 
        y(t0+i tau) = sol[  floor( n*(1+(t0 +i tau)/tau )  ]    
        on the interval [t0 + i*tau, t0 + (i+1)*tau]
     # update the global solution
     for j in range(n):
        sol[j+(i+1)*n] = lsol[j]        
# return the solution starting from t0
return sol.subarray(n, sol.length-n)

Some examples

The Delay-Action Oscillator (DAO) is a simple model of El Nino. The model is described in detail in the paper El Nino and the Delay Action Oscillator. It consists of a single delay ordinary differential equation of this form:

(5)T=TT 3αT(tδ) T' = T - T^{3} - \alpha T(t-\delta)

where TT is the deviation from the “normal” average temperature, while α\alpha and δ\delta are constants.

Below we present some numerical solutions for the DAO with α\alpha equal to 0.750.75 and different values for δ\delta.

Download code

Plots and numerical solutions presented here are generated using a homemade Java implementation of the step method. You can download the code here.

Available software

Pydelay is a Python library for solving DDE which translates DDEs into C code. Noise can be included in the simulations, too.


See also stochastic delay differential equation, and try:

Numerical methods for deterministic delay differential equations are explained here:

  • Alfredo Bellen, Marino Zennaro: Numerical methods for delay differential equations. (ZMATH)