Kuramoto model as DDE: Results don't converge when constant_lags are specified

[duodenum96]

Hello,

Here is my problem: I am trying to simulate the Kuramoto model with time delays but delays come from a matrix D where ij-th element describes the time delay between oscillator i and oscillator j. Moreover, oscillators are connected with a weighted connectivity matrix C. The model is written as

d(theta_i) / dt = omega_i + K sum_over_j( C[i, j] sin( theta_j(t - D[i, j] - theta_i(t)) )

Where omega_i is the natural frequency and thetas are phases of oscillators

Here is how I coded the function:

function kuramoto_dde!(dtheta, theta, h, p, t)
    omega, C, K, D = p

    for i in eachindex(omega)
        sum_coup = 0.
        # Do the sigma term (second term in the model)
        for j in eachindex(omega)
            if D[i, j] != 0 # Some of the connectivities are missing so delay applies only for available ones
                thetaj = h(p, t - D[i, j])[j]
            else
                thetaj = theta[j]
            end
            sum_coup += C[i, j] * sin(thetaj - theta[i])
        end
        # Write the full equation
        dtheta[i] = omega[i] + K * sum_coup
    end
    nothing
end

The code to run the simulations:

using MAT, DifferentialEquations, LinearAlgebra
# These MATLAB files contain matrices C and D
C = matread(raw"C:\Users\duodenum\Desktop\brain_stuff\pr_basics\data\chaudhuri_sc.mat")["J"]
D = matread(raw"C:\Users\duodenum\Desktop\brain_stuff\pr_basics\data\chaudhuri_D.mat")["D"]

n = size(C)[1]
dt = 0.002
tspan = (0., 300.)
timevector = 0:dt:tspan[2]

# Define the history function
h(p, t) = rand(n) * 2. * pi

# Omegas come from a gaussian distribution with mean omega_0 and variance dispersion
omega_0 = 10.
dispersion = 9.
omega = 2. * pi .* (omega_0 .+ sqrt(dispersion) .* (randn((n, 1))))

# Randomly distributed initial conditions
theta_0 = 2. * pi * rand(n)

K = 1000
p = (omega, C, K, D)

prob = DDEProblem(kuramoto_dde!, theta_0, h, tspan, p; constant_lags = D)

alg = MethodOfSteps(Rosenbrock23())
sol = solve(prob, alg; saveat = timevector)
println("SOLUTION COMPLETE")

Which gives the following error:

Warning: dt(0.0) <= dtmin(2.220446049250313e-16) at t=0.0, and step error estimate = 1.0. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase C:\Users\duodenum\.julia\packages\SciMLBase\GmToj\src\integrator_interface.jl:599

On the other hand, if I don't specify constant lags:

prob = DDEProblem(kuramoto_dde!, theta_0, h, tspan, p)

alg = MethodOfSteps(Rosenbrock23())
sol = solve(prob, alg; saveat = timevector)
println("SOLUTION COMPLETE")

I get the solution but it is very slow. I tried using constant_lags = D[:] or constant_lags = D[:]' as well but none of them worked. Is there a way to specify constant_lags in this situation without messing up the results?

Thanks for the great work with this package (and the sciML ecosystem in general!) Yasir

[ChrisRackauckas]

As described at https://stackoverflow.com/questions/77195203/kuramoto-model-as-dde-results-dont-converge-when-constant-lags-are-specified, the model is ill-posed so it needs to be fixed to have an actually defined history.

[duodenum96]

Thanks for the solution. Just writing the following in case anyone has the same problem: turns out the delay matrix has 0s inside it and the constant_lags option on differentialequations doesn't like 0s on constant_lags. The following modification solved the problem:

lags = D_over_v[:]
deleteat!(lags, lags .== 0)
prob = DDEProblem(kuramoto_dde!, theta_0, h, tspan, p; constant_lags = lags)