{{ message }}

Instantly share code, notes, and snippets.

# Zulko/ddeint.py

Created Oct 22, 2013
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 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 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 commented Dec 3, 2015

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

### petehellyer commented Mar 7, 2016

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

### 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 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.