I'm using the JiTCDDE module to solve delayed differential equations. My work requires me to let the model evolve for a while and then perturb it. To that end I'm trying to use jitcdd.purge_past()
followed by jitcdde.add_past_points()
. Issue is, the latter requires me to provide not only the values at the times I input, but also the values of the derivative.
Is there a way to get those from the jitcdde
instance itself, or do I need to manually and clunckily calculate them?
EDIT: more info
My system consists of two non linear oscillators coupled together (they are not phase oscillators). I let the system evolve for a while until it reaches a steady state and then perturb it by shifting one of the two oscillators in time a bit. The amount I shift it is calculated as a percentage of the period of oscillation, what amounts to an effective phase shift.
My system is 12-dimensional and I'm struggling to get a minimal working example going, so here's a non-working mockup of the code:
f = [...]
DDE = jitcdde(f)
DDE.constant_past([1]*len(f))
DDE.step_on_discontinuities()
initial_integration_time = 30
DDE.integrate(initial_integration_time)
After doing this, I'd like to perform the perturbation that should be, lets say 10% of a period, assume I know the stationary period length to be T
. What I'm trying to do is to use get_state
to get the current state of the system and derivative, and the state and derivative perturbation
time units ago. Then, I can construct a new set of anchor
s that I can pass to DDE.add_past_points
to start integrating from there.
T = ...
perturbation = 0.1*T #in units of time
state = DDE.get_state()
# get position and derivative of first oscilator
# since the system is 12-D, the first six correspond to osc1
osc1_x = [state[1][6:] for s in state]
osc1_dx = [state[2][6:] for s in state]
# here, do the same for osc2 at DDE.t - perturbation
To answer the question as posted, to get the derivative from a jitcdde instance all you need to do is call the get_state
method. Assuming DDE
is your jitcdde
instance which you already integrated:
state = DDE.get_state()
derivatives = [s[2] for s in state]
This works because get_state
will return a Cubic Hermite Spline instance that behaves basically like a list of tuples (plus some cool methods) where each tuple contains a time, the state of the system at that time and the derivatives of the system at that time, see this entry in the docs for more info.
To solve my problem specifically, if any passerby happens to care, I did the following (code in the same mock-up style as the question):
# define and initialize system
f = [...]
DDE = jitcdde(f)
DDE.constant_past([1]*len(f))
DDE.step_on_discontinuities()
initial_integration_time = 30
T = ... # period length
perturbation = 0.1*T # 10% of period, in units of time
# get state of system at perturbation time units before the
# end of the initial integration
DDE.integrate(initial_integration_time-perturbation)
state1 = DDE.get_state.copy()
# get state after initial integration is done
DDE.integrate(initial_integration_time)
state2 = DDE.get_state()
# generate new anchors from these set of states
eps = 1e-6
anchors = []
for i, s in enumerate(state2[::-1]): #iterate in reverse
# stop if I reached the end, there's no telling
# which states is gonna be longer
if i >= len(state1):
break
# calculate new anchor at the time of the anchor I'm looking at
# shifted by the perturbation length
tt = s[0] - perturbation_duration
state1.interpolate_anchor(tt)
_, x1, dx1 = state1[state1.last_index_before(tt+eps)]
x_new = np.hstack( (x1[:6], s[1][6:]) )
dx_new = np.hstack( (dx1[:6], s[2][6:]) )
anchors.append( (tt, x_new, dx_new) )
# add this new past to the integrator and evolve the system
DDE.purge_past()
DDE.add_past_points(anchors)
DDE.step_on_discontinuities()
# now I can evolve the system as I please
This code has some caveats, like ensuring the perturbation is not too big and that the states are long enough, but this is the basic idea behind my solution. I'm not gonna explain what it does in detail, since I don't think it'd be instructive to anyone.