Question

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.

Was it helpful?

OTHER TIPS

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.

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.

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top