Getting derivative value from JiTCDDE instance

Question

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 anchors 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

Answer 1

The question

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.

---

My problem

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.