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)]) |
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
Same error here.
Did anyone solve this issue?
Best and Thx,
Felix