# Euler integration of delay differential equation with circular buffer

I'm seeking a way to streamline the updating buffer process. The code provided addresses a coupled differential equation with a delay. However, it's worth noting that it doesn't employ a circular buffer as intended, instead opting for a rolling buffer approach which is quite time-consuming.

Therefore, the primary objective is to devise a memory-efficient solution for solving a system of delay differential equations without necessitating the storage of the entire time series.

\begin{align} dx &= -x(t)y(t-\tau_0) + y(t-\tau_1) \\ dy &= x(t)y(t-\tau_0) - y(t) \\ dz &= y(t) - y(t-\tau_1) \end{align}

Refer to the Euler method.

import matplotlib.pyplot as plt
import numpy as np
np.random.seed(2)

def fn(y, D):
f0 = -y[0, -1] * y[1, -D[0]] + y[1, -D[1]]
f1 = y[0, -1] * y[1, -D[0]] - y[1, -1]
f2 = y[1, -1] - y[1, -D[1]]
return np.array([f0, f1, f2])

def intg():

dt = 0.01
y0 = np.array([5.0, 0.1, 1.0])
L = np.array([1.0, 10.0])
t = np.arange(0.0, 100, dt)
nstep = len(t)
tcut = 5
decimate = 10
nbuffer = int((t[-1]-tcut)//(dt*decimate) + 1)
D = np.round(L/dt).astype(np.int64)
md = np.max(D)                  # maximum delay
y = np.tile(y0, (md, 1)).T      # store buffer
Y = np.zeros((len(y0), nbuffer))
store_count = 0

for it in range(nstep):
curr_t = it * dt
y_new = fn(y, D) * dt + y[:, -1]

# update buffer  (Time consuming part)
y[:, :-1] = y[:, 1:]
y[:, -1] = np.reshape(y_new, (3))

if curr_t >= (tcut):
if (it % decimate) == 0:
Y[:, store_count] = np.reshape(y_new, (3))
store_count += 1
return Y, t[t >= tcut][::decimate]

Y, t = intg()

print(Y.shape, t.shape)
plt.plot(Y.T)
plt.show()


buffer_index = (buffer_index + 1) % md
y[:, buffer_index] = y_new


and correspondingly changing the equations like this: y[0, -D[0]) -> y[0, (t-D[0]) % md)

• Could you spell out the DE you are integrating? This exchange supports LaTeX.
– vnp
Feb 22 at 16:21
• @J_H The main idea is to use a memory efficient method to solve a delay differential equation, without needing to store the whole time series. The given system of equations is a simple toy example, These kind of equations can be used in 1. Population dynamics 2. Disease outbreaks 3. Traffic flow and 4.Neural networks Feb 29 at 10:19

Your integration method is not reusable; you need to make it reusable. The parameters are fairly easy to add and hint. For the integrand callable you can use a Protocol.

update buffer (time-consuming part)

This is possible but not definite. You haven't shown any evidence of profiling, and you need to profile; or else any explanation of a performance mechanism will be a wild guess. This will only be useful if you also break down your function into sub-routines.

You label this a circular buffer, but it isn't a circular buffer - you completely overwrite the buffer on every iteration. For it to be a proper circular buffer, this:

        y[:, :-1] = y[:, 1:]
y[:, -1] = y_new


needs to be rewritten so that

• only one column is written,
• the column write index increments on each iteration, and
• the column write index and the D indices are the subject of a modulus.

Stop using the implicit interface to matplotlib and use the object interface instead. Also, you need to be plotting against your time series instead of making matplotlib infer indices.

Covering some of the above, without the modular re-write, we have

from typing import Protocol

import matplotlib.pyplot as plt
import numpy as np

def my_integrand(y: np.ndarray, D: np.ndarray) -> np.ndarray:
f0 = -y[0, -1] * y[1, -D[0]] + y[1, -D[1]]
f1 = y[0, -1] * y[1, -D[0]] - y[1, -1]
f2 = y[1, -1] - y[1, -D[1]]
return np.array((f0, f1, f2))

class Integrand(Protocol):
def __call__(self, y: np.ndarray, D: np.ndarray) -> np.ndarray:
...

def euler_integrate(
fn: Integrand,
y0: np.ndarray,
L: np.ndarray,
tcut: float,
t1: float,
t0: float = 0,
dt: float = 0.01,
decimate: float = 10,
) -> tuple[
np.ndarray,  # Y
np.ndarray,  # t
]:
t = np.arange(t0, t1, dt)
nstep = len(t)
nbuffer = int((t[-1] - tcut)//(dt * decimate)) + 1
D = np.round(L / dt).astype(np.int64)
md = np.max(D)  # maximum delay
y = np.tile(y0, (md, 1)).T  # store buffer
Y = np.zeros((len(y0), nbuffer))
store_count = 0

for it in range(nstep):
curr_t = it*dt
y_new = fn(y, D)*dt + y[:, -1]

# update buffer (time-consuming part)
y[:, :-1] = y[:, 1:]
y[:, -1] = y_new

if curr_t >= tcut and it % decimate == 0:
Y[:, store_count] = y_new
store_count += 1

return Y, t[t >= tcut][::decimate]

def demo() -> None:
Y, t = euler_integrate(
fn=my_integrand,
y0=np.array((5, 0.1, 1)),
L=np.array((1., 10.)),
t1=100,
tcut=5,
)
print(Y.shape, t.shape)

fig, ax = plt.subplots()
for series in Y:
ax.plot(t, series)
plt.show()

if __name__ == '__main__':
demo()


• good points, probably here was not a proper place to ask, but I was mainly looking for resolving the expensive part of copying the whole buffer and replace with a circular buffer. Feb 29 at 10:07

## Naming

Some of your variables have meaningful names, but others are not too descriptive. It is fine to have some variables with short names, but in those cases, you should add comments to describe them. For example, dt, Y, t, etc.

You could also use more descriptive names for your functions: fd, etc.

## Output

You could print out the description of the text output. I see this when I run the code:

(3, 950) (950,)


It would be helpful to see what those numbers refer to.

You could add descriptive labels labels to the plot:

plt.plot(Y.T)
plt.title("stuff")
plt.xlabel("x stuff")
plt.ylabel("y stuff")
plt.show()


## Unused code

It seems the following code is not needed:

np.random.seed(2)


It is a good idea to remove unused code.