Create a gist now

Instantly share code, notes, and snippets.

A simple Delay Differential Equation solver written in Python, using the solving capabilities of the Scipy package.
# 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)])
@weiszr
weiszr commented Dec 9, 2014

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

@fkemeth
fkemeth commented Dec 3, 2015

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

@petehellyer

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

@julianeagu
julianeagu commented Oct 31, 2016 edited

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.

@Ziaeemehr
Ziaeemehr commented Nov 15, 2016 edited

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment