__del__( self )

Eaten by the Python.

Delay Differential Equations in Python

| Comments

I wrote ddeint, a simple module/function for solving Delay Differential Equations (DDEs) in Python. It is not very fast, but very flexible, and coded in just a few lines on top of Scipy’s differential equations solver, odeint.

Say you have a delay differential equation like this:

where $F(y, t)$ can involve delayed values of $y$, of the form $y(t-d)$.

To solve this DDE system at points t=[t1, t2 ...] you would just write

1

y = ddeint(F, g, t) # returns [y(t1), y(t2) ...]

A simple example

Let us start with a DDE whose exact solution is known (it is the sine function), just to check that the algorithm works as expected:

Here is how we solve it with ddeint:

1
2
3
4
5
6
7
8
9
10
11
12
13

from pylab import *
from ddeint import ddeint

model = lambda Y,t : Y(t - 3*pi/2) # Model
tt = linspace(0,50,10000) # Time start, time end, nb of points/steps
g=sin # Expression of Y(t) before the integration interval
yy = ddeint(model,g,tt) # Solving

# PLOTTING
plot(tt,yy,c='r',label="$y(t)$")
plot(tt,sin(tt),c='b',label="$sin(t)$")
set_ylim(ymax=2) # make room for the legend
legend()

The resulting plot compares our solution (red) with the exact solution (blue). See how our result eventually detaches itself from the actual solution as a consequence of many successive approximations ? As DDEs tend to create chaotic behaviors, you can expect the error to explode very fast. As I am no DDE expert, I would recommend checking for convergence in all cases, i.e. increasing the time resolution and see how it affects the result. Keep in mind that the past values of Y(t) are computed by interpolating the values of Y found at the previous integration points, so the more points you ask for, the more precise your result.

An example with parameters

You can set the parameters of your model at integration time, like in Scipy’s ODE and odeint. As an example, imagine a chemical product with degradation rate $r$, and whose production rate is negatively linked to the quantity of this same product at the time $(t-d)$:

We have three parameters that we can choose freely. For $K = 0.1$, $d = 5$, $r = 1$, we obtain oscillations !

1
2
3
4
5
6
7
8

# MODEL, WITH UNKNOWN PARAMETERS
model = lambda Y,t, k,d,r :  1/(1+(Y(t-d)/k)**2) - r*Y(t)
g = lambda t:0 # history before t=0

tt = linspace(0,50,10000)
yy = ddeint(model, g, tt, fargs=(0.1, 5, 1)) # K=0.1, d=5, r=1

plot(tt,yy,lw=2)

Example with several variables

The variable Y can be a vector, which means that you can solve DDE systems of several variables. Here is a version of the famous Lotka-Volterra two-variables system, where we introduce a delay $d$. For $d=0$ the system is a classical Lotka-Volterra system ; for $d\neq 0$ the system undergoes an important amplification:

1
2
3
4
5
6
7
8
9
10
11
12
13
14

def model(Y,t,d):
    x,y = Y(t)
    xd,yd = Y(t-d)
    return array([0.5*x*(1-yd), -0.5*y*(1-xd)])

g = lambda t : array([1,2])
tt = linspace(2,30,20000)

for d in [0, 0.2]:
    yy = ddeint(model,g,tt,fargs=(d,))
    # WE PLOT X AGAINST Y
    plot(yy[:,0], yy[:,1], lw=2, label='delay = %.01f'%d)

legend() # display the legend

Example with a non-constant delay

In this last example the delay depends on the value of $y(t)$ :

1
2
3
4
5

model = lambda Y,t:  -Y( t-3*cos( Y(t) )**2 )
tt = linspace(0, 30, 2000)
yy = ddeint(model, lambda t:1, tt)

plot(tt, yy)

Comments