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.

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 \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-\tau_{1}), \ldots , y(t-\tau_{m}) )$

(2)$y(t) = g(t), \; t_{0} - \tau \lt t \leq t_{0}$

where $g$ is given and represents the values of $y$ for $t_{0} - \tau \lt t \leq t_{0}$. More about the theory of delay differential equations is discussed here.

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-\tau))$

(4)$y(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 \tau]$ with $k \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}+\tau]$, $[t_{0} + \tau, t_{0}+2 \tau]$ and so on. In fact, when we solve the problem for $[t_{0} + i \tau, t_{0}+(i+1) \tau]$, we already know the values of $y(t)$ for $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)

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' = T - T^{3} - \alpha T(t-\delta)$

where $T$ 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.75$ and different values for $\delta$.

- First, the solution for $\alpha=0.75, \delta=1.0$ (download data here):

- Second, the solution for $\alpha=0.75, \delta=2.0$ (download data here):

- Third, the solution for $\alpha=0.75, \delta=4.0$ (download data here):

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

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:

- Bruce E. Shapiro, Topics in numerical analysis, course notes, 2007

Numerical methods for deterministic delay differential equations are explained here:

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