A simple Delay Differential Equation solver written in Python, using the solving capabilities of the Scipy package.

ddeint.py

# REQUIRES PACKAGES Numpy AND Scipy INSTALLED
import numpy as np
import scipy.integrate
import scipy.interpolate

class ddeVar:
    """ special function-like variables for the integration of DDEs """

    def __init__(self,g,tc=0):
        """ g(t) = expression of Y(t) for t<tc """

        self.g = g
        self.tc= tc
        # We must fill the interpolator with 2 points minimum
        self.itpr = scipy.interpolate.interp1d(
            np.array([tc-1,tc]), # X
            np.array([self.g(tc),self.g(tc)]).T, # Y
            kind='linear', bounds_error=False,
            fill_value = self.g(tc))

    def update(self,t,Y):
        """ Add one new (ti,yi) to the interpolator """

        self.itpr.x = np.hstack([self.itpr.x, [t]])
        Y2 = Y if (Y.size==1) else np.array([Y]).T
        self.itpr.y = np.hstack([self.itpr.y, Y2])
        self.itpr.fill_value = Y

    def __call__(self,t=0):
        """ Y(t) will return the instance's value at time t """

        return (self.g(t) if (t<=self.tc) else self.itpr(t))

class dde(scipy.integrate.ode):
    """ Overwrites a few functions of scipy.integrate.ode"""

    def __init__(self,f,jac=None):

        def f2(t,y,args):
            return f(self.Y,t,*args)
        scipy.integrate.ode.__init__(self,f2,jac)
        self.set_f_params(None)

    def integrate(self, t, step=0, relax=0):

        scipy.integrate.ode.integrate(self,t,step,relax)
        self.Y.update(self.t,self.y)
        return self.y

    def set_initial_value(self,Y):

        self.Y = Y #!!! Y will be modified during integration
        scipy.integrate.ode.set_initial_value(self, Y(Y.tc), Y.tc)

def ddeint(func,g,tt,fargs=None):
    """
    Similar to scipy.integrate.odeint. Solves a Delay differential
    Equation system (DDE) defined by ``func`` with history function ``g``
    and potential additional arguments for the model, ``fargs``.
    Returns the values of the solution at the times given by the array ``tt``.
    
    Example:
    --------
    
    We will solve the delayed Lotka-Volterra system defined as
    
    For t < 0:
    x(t) = 1+t
    y(t) = 2-t
    
    For t > 0:
    dx/dt =  0.5* ( 1- y(t-d) )
    dy/dt = -0.5* ( 1- x(t-d) )
    
    Note that here the delay ``d`` is a tunable parameter of the model.
    
    ---

    import numpy as np

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

    g = lambda t : np.array([1+t,2-t]) # 'history' at t<0
    tt = np.linspace(0,30,20000) # times for integration
    d = 0.5 # set parameter d 
    yy = ddeint(model,g,tt,fargs=(d,)) # solve the DDE !
     
    """

    dde_ = dde(func)
    dde_.set_initial_value(ddeVar(g,tt[0]))
    dde_.set_f_params(fargs if fargs else [])
    return np.array([g(tt[0])]+[dde_.integrate(dde_.t + dt)
                                 for dt in np.diff(tt)])

Hello,
I have been trying to play with your solver, but I get the following error: index 2 is out of bounds for axis 0 with size 2. The full error is copied below.

Do you have any idea what the problem is?

Cheers,

Robert

Full error:
Traceback (most recent call last):
File "text.py", line 13, in
yy = ddeint(model,g,tt,fargs=( 0.1 , 5 , 1 )) # K = 0.1, d = 5, r = 1
File "/Users/weiszr/Dropbox/apps/Computable/ddeint.py", line 97, in ddeint
for dt in np.diff(tt)])
File "/Users/weiszr/Dropbox/apps/Computable/ddeint.py", line 46, in integrate
scipy.integrate.ode.integrate(self,t,step,relax)
File "/Users/weiszr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/integrate/_ode.py", line 388, in integrate
self.f_params, self.jac_params)
File "/Users/weiszr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/integrate/_ode.py", line 737, in run
args[5:]))
File "/Users/weiszr/Dropbox/apps/Computable/ddeint.py", line 40, in f2
return f(self.Y,t,_args)
File "text.py", line 6, in
model = lambda Y,t,k,d,r : 1/(1+(Y(t-d)/k)__2) - r_Y(t)
File "/Users/weiszr/Dropbox/apps/Computable/ddeint.py", line 32, in call
return (self.g(t) if (t<=self.tc) else self.itpr(t))
File "/Users/weiszr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/interpolate/polyint.py", line 79, in call
y = self._evaluate(x)
File "/Users/weiszr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/interpolate/interpolate.py", line 478, in _evaluate
y_new = self._call(self, x_new)
File "/Users/weiszr/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/interpolate/interpolate.py", line 441, in _call_linear
y_hi = self._y[hi]
IndexError: index 2 is out of bounds for axis 0 with size 2

Same error here.
Did anyone solve this issue?
Best and Thx,
Felix

I also have the same problem, has anyone been able to solve?

Following the solution listed here and here, one needs to add

self.itpr._y = self.itpr._reshape_yi(self.itpr.y)

after line 28. This solved the issue for me.

I compared the result of ddeint and pydelay and here is the code:

""" Reproduces the sine function using a DDE """
from pylab import *
from ddeint import ddeint
from pydelay import dde23
from timeit import default_timer as timer
from numpy import *

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

end = timer()
print "elapsed time for ddeint : %.4f seconds" % (end - start)

plot(tt,yy,'r',label="$ddeint$")
plot(tt,sin(tt),'b',label="$sin(t)$")
ylim([-1.5,2])

start = timer()
par = {'tau': 1.5*np.pi}
eqns = {'x': 'sin(t-tau)'}
dde = dde23(eqns = eqns, params = par )
dde.set_sim_params(tfinal=50, dtmax=0.1)
histfunc = {
    'y': lambda t : sin(t)
    }
dde.run()
sol = dde.sample(0,50,2e-3)
y1 = sol['x']
t1 = sol['t']
end = timer()
print "elapsed time for pydelay : %.4f seconds" % (end - start)

plot(t1,y1,'k-',label='pydelay')

legend()
show()

And here the results:

elapsed time for ddeint : 6.5165 seconds
elapsed time for pydelay : 0.0186 seconds

I should also mention that the resulting curve from pydelay was very close to the analytic solution and ddeint digressed.