Created
October 22, 2013 06:24
-
-
Save Zulko/7096011 to your computer and use it in GitHub Desktop.
A simple Delay Differential Equation solver written in Python, using the solving capabilities of the Scipy package.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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)]) |
Same error here.
Did anyone solve this issue?
Best and Thx,
Felix
I also have the same problem, has anyone been able to solve?
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
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