# The Azimuth Project Numerical methods for delay differential equations

## Idea

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 \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.

## 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-\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)


### 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' = 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.

## Available software

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