CODEWITHSUNDEEP
juliadifferential-equationsdifferentialequations.jl
duodenum
duodenum

Reputation: 95

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

Copy-pasting from the issue I opened on github (https://github.com/SciML/DifferentialEquations.jl/issues/983).

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

Upvotes: 0

Views: 182

Answers (2)

duodenum
duodenum

Reputation: 95

Thanks for various tips. Chris's answer makes perfect sense. Posting this just in case somebody else has a similar problem: In my case, the delay matrix D had 0s in the diagonal and it turns out the lags vector doesn't like 0s inside it. So the following modification solved the error:

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

Upvotes: 0

Chris Rackauckas
Chris Rackauckas

Reputation: 19132

This is ill-posed. h(p, t) = rand(n) * 2. * pi gives a different random number each time, since it's just a function calling randn. The history is is almost surely not continuous or differentiable and isn't even well-defined (it does not evaluate to the same value each time). The DDE solution is not defined in any sensible mathematical way. The adaptivity is probably trying really hard to figure out why it cannot get a single unique solution and gives up because there is no solution to the model as posed.

Upvotes: 2

Related Questions