Nice.

I usually work on a system of DDE with constant and heterogeneous delays.

Some thing like this for a system of Wilson-Cowan equations:

or

delayed couple Kuramoto oscillators that I already worked on.

For a toy example letâ€™s consider the following DDE and step by step we can add more complexity to our solver.

```
dy0dt = -y0 + y1(t-1) + y1(t-10)
dy1dt = y0 + y1(t-1) - y1
dy2dt = y1 - y1(t-10)
```

we suppose the delays are integer multiple of time step `dt`

(we can release this restriction later).

For numerical integration we evaluate the state of the system at discrete time grids with equal space (for fixed time step) or unequal ones (adaptive methods).

If delays are integer multiple of time step (the delays need to be larger than the time step)

we can have access to previous states of the system only by using the proper index like `y[n-delay_index0][1]`

at the following equations.

For **Euler** scheme :

There is no need to touch the standard scheme.

The integrator look like this:

```
def euler(f, a, b, h):
'''
f: (function) RHS of the DDE
a : [float] start time
b : [float] end time
h : [float] time step
'''
n_iteration = (int)((b-a)/h)
n = len(y)
for i in range(nstart, nstart+n_iteration):
tmp = h * f(y[i], 0.0, i)
y.append(tmp + y[i])
t.append((1+i)*dt-maxdelay)
```

The equations are :

```
def dde(x, t, n):
# define dde system of equations
f0p = -y[n][0] * y[n-delay_index0][1] + y[n-delay_index1][1]
f1p = y[n][0] * y[n-delay_index0][1] - y[n][1]
f2p = y[n][1] - y[n-delay_index1][1]
return np.asarray([f0p, f1p, f2p])
```

where

```
delays = [1.0, 10.0]
delay_index0 = int(delays[0]/dt)
delay_index1 = int(delays[1]/dt)
maxdelay = max(delays)
```

we also need to prepare some history for the beginning of the simulation.

at least `int(maxdelay/dt)`

steps. The history values affect the numerical value of the next steps. Usually we know the history or just put a random or fixed value for that.

```
def set_history(hist, nstart, maxdelay):
N = len(hist)
for i in range(nstart+1):
t.append(-(nstart - i) / float(nstart) * maxdelay)
# x: nstart x N
for i in range(nstart+1):
y.append(hist)
```

this example use a list of list to record state of the system at each time step which is slow and can be replaced with 2d numpy array.

To make program more memory efficient we probably can store only `int(maxdelay/dt)`

time step.

If we want to have arbitrary delays, then we donâ€™t have the previous points necessarily on the evaluated points on the time grid. So we need to interpolate the state of the system on those points. `np.interp`

is a choice.

The post will be too long to cover other integration schemes. May be Itâ€™s better to send you some test and benchmarks. I will be available if there be any ambiguity.

This is the result:

```
dt = 0.001
y0 = [5.0, 0.1, 1.0]
delays = [1.0, 10.0]
delay_index0 = int(delays[0]/dt)
delay_index1 = int(delays[1]/dt)
maxdelay = max(delays)
nstart = delay_index1
t_sim = 100
y = []
t = []
set_history(y0, nstart, maxdelay)
euler(dde, 0, t_sim, dt)
```

The example comes from

dde23 tutorial for matlab page 16. If you like to check the results.

Reference: look the book appendix

Delay-Coupled Complex Systems and Applications to Lasers

Flunkert , Valentin

I agree with any of suggested syntax.

cheers,

Abolfazl