3

Basically... I need a way to include a phase shift in my differential equations. That is, I have in the definition of my system function which returns dY/dt something like Y(t-3). Like this differential equation:

dY/dt = a*Y(t) + b*Y(t-tau)

Now if I try to write this as the system definition function for passing to scipy.odeint, I am lost:

def eqtnSystem(A,t):
    Y   = A
    a   = 1
    b   = 5
    tau = 3
    return a*Y + b*???       # how do I Y(t-tau) ?

That's basically it. I really hope there is an easy answer, but I couldn't seem to track one down.

Specifically... I am attempting to numerically calculate the solution for the system defined by the following function:

def etaFunc(A,t): 
    #...definition of all those constants is here...
    return array([(gamma[0,0]*xi(t-theta[0])[0] - eta[0] + zeta[0])/tau[0],\
           (gamma[1,1]*xi(t-theta[1])[1] - eta[1] + zeta[1])/tau[1],\
           (gamma[2,2]*xi(t-theta[2])[2] - eta[2] + zeta[2])/tau[2],\
           (   beta[3,0]*pastEta(t-theta[3])[0] \
             + beta[3,1]*pastEta(t-theta[4])[1] \
             + beta[3,2]*pastEta(t-theta[5])[2] -eta[3]+ zeta[3])/tau[3],\
           (   beta[4,3]*pastEta(t-theta[6])[3] \
             + beta[4,2]*pastEta(t-theta[7])[2] - eta[4] + zeta[4])/tau[4]])

This function is then later given to odeint like this:

ETA = integrate.odeint(etaFunc,initCond,time)

and then I can get out each individual component of ETA (like eta_0) like this: ETA[:,0].

The problem I am having here, is with pastEta(t-theta[?]). For right now, this is a function which attempts to find already calculated values of eta (for when start_time < t-theta[?] < t and theta[?] > 0. This isn't working very well.

I see in this case I could find each component of eta individually and then get calculated values for previously calculated eta components (eta_0,eta_1,eta_2) to calculate eta_3 and similarly for eta_4, but this is not ideal since it takes away the ability for me to 'plug-and-play' any general formulas.

3 Answers 3

3

There are a number of existing libraries and examples for doing this.

http://www.google.fi/search?q=python+delay+differential+equation gives me:

2
  • Thanks so much! I missed all these somehow; I'll call it 'delay' from now on.
    – 7yl4r
    Apr 4, 2013 at 19:23
  • +1: These methods appear to use the approach that I described in my answer, but save the need to build it oneself, so they are much to be preferred.
    – Simon
    Apr 4, 2013 at 20:01
2

Delays aren't exactly linear functions. The usual step delay is represented in Laplace domain as e**(a*s)/s, where a is the delay.

What this means is that "normal" ODE solvers won't work unless you have some workaround. Usually this workaround isn't very easy to do, since for stiff problems you usually can't interpolate with a good enough approximation.

Anyway, one of the solutions is using the libraries posted in the other answers.

Another solution is doing it symmbolically (if you can, you might try SymPy).

A third solution is storing the past results and interpolating to find the exact past you need (might not be good enough).

A fourth solution might be what is recommended by some simulators docs: use c2d() and simulate the whole model in discrete time and store past variables in an list/array (no interpolation, but you might need to use small steps for better accuracy).

A fifth solution is using Padé approximation to represent your model's delay (might working depending on your case). There's a pade() function in python-control to approximate exactly this.

1

One way to do this with integrate.odeint() would be to run integrate.odeint() for many short time intervals between your starting time and your ending time, storing the time value and the output Y value after each short interval in lists. That would let you interpolate the Y value in the lists using scipy.interpolate.interp1d(), for instance, each time you needed Y(t-3).

You only end up with an approximate value for Y(t-3) if you do it this way, of course, but if the time intervals are close enough together, this approach might be satisfactory for you. After all, the Y(t) values calculated by numerical ODE solvers are approximate too.

4
  • Ah yes... this could work so long as the phase shift is negative. I'm surprised odeint doesn't have any built-in functionality for this though. It seems like it would be a fairly common problem.
    – 7yl4r
    Apr 4, 2013 at 15:11
  • @7yl4r: The problem is harder if it involves a positive phase shift because the ODE is no longer solvable as an initial value problem and different methods are required.
    – Simon
    Apr 4, 2013 at 19:55
  • @7yl4r: You could now accept your preferred answer by checking or ticking the one that you've found most useful. If you don't have a strong preference for either answer, I'd suggest accepting pv's answer because it is more immediately useful than mine.
    – Simon
    Apr 6, 2013 at 20:23
  • Thanks Simon. I found them both very useful - yours for the concept, and pv's for the implementation.
    – 7yl4r
    Apr 10, 2013 at 2:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.