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

Using `get_state` should allow for most perturbations of the past. Can you give us some detail why it does not work for you? Also, what kind of perturbation do you wish to perform, e.g., a pulse, a step, etc.? — Wrzlprmft, Nov 04 '21 at 13:30
It's not that `get_state` didn't work, it's just that I didn't look deep enough into it to understand how to use it. I'll take a look at it and see if I can figure it out. The issue is not that I don't know how to get the past, is that I don't know hot to get f evaluated in the past, aka, the past derivatives. My system is two very non liner coupled oscillators and the perturbation is a phase shift of one of them. — Puff, Nov 04 '21 at 17:13
Please [edit] the information on your perturbation into your question. If you additionally specify whether you have phase oscillators (i.e., oscillators where the phase is a dynamical variable) or not, I can answer your question. — Wrzlprmft, Nov 04 '21 at 20:12
I didn't add a specific example because I'm struggling to come up with a minimal working one, and I think that adding my 12-dimensional system here is not gonna be that helpful. I'll try to write something up, but meanwhile, I'll make the question more specific. — Puff, Nov 04 '21 at 20:16

score 1 · Accepted Answer · edited Dec 13 '21 at 15:42

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.

You can probably do the perturbation a bit more elegantly by creating separate splines for each oscillator (with one of them time shifted) and then [joining](https://chspy.readthedocs.io/en/latest/#_chspy.join) them and finally [truncating](https://chspy.readthedocs.io/en/latest/#_chspy.CubicHermiteSpline.truncate) the result. — Wrzlprmft, Nov 05 '21 at 18:51
Ahh, I see! If I implement it, I'll update the answer. Thanks! — Puff, Nov 08 '21 at 01:02
How do you suggest I split the coordinates that correspond to each oscillator so I can shift and join them? I don't see a more elegant way than what I did, which is not elegant at all, mind you. — Puff, Nov 08 '21 at 14:37

Getting derivative value from JiTCDDE instance

1 Answers1

The question

My problem