Extracting derivatives for delay differential equations with JiTCDDE

Question

I'm simulating the time-delayed Kuramoto model using JiTCDDE and I am having some issues using get_state() to extract the derivatives through the time evolution.

I'm modelling the time-delayed Kumamoto model on a digraph encoded by an adjacency matrix A:

def kuramotos():
    sols = np.zeros(n)
    for i in range(n):
         yield omegas[i] + k * sum(
             sin(y(j, t - τ) - y(i))
            for j in range(n) if A[j,i])

I'm using JiTCDDE to simulate the time evolution according to:

I = jitcdde(kuramotos, n=n)

I.set_integration_parameters(rtol=0, atol=1e-5)
I.constant_past(random.uniform(0, pi, n), time=0)
I.integrate_blindly(max(τ), 0.001)
ts = np.linspace(0,100,1001)

p = [] #phases
dp = [] #phase velocities

for time in I.t+np.arange(0,100.1,0.1):
    
    p.append(list(I.integrate(time) % (2 * pi)))
    state = I.get_state() # returns time, phases, frequency
    dpi = [s[2] for s in state]
    dp.append(dpi)

From what I gather, [s[2] for s in state] should be returning the values of the derivatives at each time. So I would expect to have a set of n phase velocities at each time step. In reality, what I'm getting is arrays of varying shapes at each time. For example, for n = 50, I would ideally be expecting a simple list of 50 values for the derivative at each timestep, but what I get is:

for i in range(0,10):
    print(np.shape(np.array(dp[i])))

Output:

(21, 50)
(7, 50)
(5, 50)
(4, 50)
(3, 50)
(3, 50)
(3, 50)
(3, 50)
(3, 50)
(3, 50)

I'm curious why the shapes of the arrays are like this, and why it seems to settle towards (3,50)? Ideally I would like to simply be getting a list of 50 values for the derivative at each step in time.

Welcome to Stack Overflow. While I could answer your question as it is, please bear in mind that it’s usually easier to post a complete non-working example, so it’s easier for others to reproduce your problem and possibly find an answer. — Wrzlprmft, Aug 21 '23 at 12:03

Wrzlprmft · Accepted Answer · 2023-08-21T12:07:39.320

get_state gives you the complete state of the integrator, i.e., what you would need to resume the integration without keeping the integrator object (I). Since DDEs depend on the past, this includes the past, i.e., multiple anchor points (see below). And since JiTCDDE is an adaptive integrator, the numbers of points used for storing the past may differ over time – which explains why the shape varies and seems to converge.

The output of get_state are the anchor points of a cubic Hermite spline, which is what JiTCDDE uses to store the past. These anchor points consist of a time, a state, and a derivative. If you want to obtain the most recent derivative, you need to extract the derivative of the last anchor point, which you can do with:

I.get_state()[-1].diff

Mind that this anchor point will not necessarily be at whatever time point you integrated to (but a bit in the future). Thus for matching results, you should extract the time and state from get_state() too. So, all in all, you want to do something like:

p = [] # phases
dp = [] # phase velocities
ts = [] #actual times

for time in I.t+np.arange(0,100.1,0.1):
    I.integrate(time)
    last_anchor = I.get_state()[-1]
    p.append( last_anchor.state % (2*pi) )
    dp.append( last_anchor.diff )
    ts.append( last_anchor.time )

Extracting derivatives for delay differential equations with JiTCDDE

1 Answers1