Kuramoto-like coupled oscillators in DifferentialEquations.jl

General Usage
,
1

Hi
I’d like to know how to implement a Kuramoto-like (coupled nonlinear oscillators) system in DifferentialEquations.jl. If you are familiar with this stuff, I am trying to produce a system proposed by Matthew,Mirollo and Strogatz Here is the paper

To my understanding, this is just a system of ODE. I produced a minimal example of it here with three oscillators(below). The thing is that they use 800 oscillators.

How one is supposed to implement (elegantly) a large system of DiffEqs in DifferentialEquations.jl ??? With Metaprogramming? Another magic?

My MWE:

Each oscillator z_{j} is just
z_{j}=(1-|z_{j}|^2+iw_{j})z_{j}+K(mean(all \,\, z_{j}) -z_{j})

function parametrized_strogatz(du,u,p,t)
#Coupling constant
K=p.K
##order parameter
du[1]=(du[2]+du[3]+du[4])/3.0   
#oscillators
du[2] = (1.0-abs2(u[2])+p.w2*im)u[2]+K*(u[1]-u[2])
du[3] = (1.0-abs2(u[3])+p.w3*im)u[3]+K*(u[1]-u[3])
du[4] = (1.0-abs2(u[4])+p.w4*im)u[4]+K*(u[1]-u[4])
end
u0 = [1.0+0.0im,1.0+0.0im,1.0+0.0im] ##starting points
pushfirst!(u0,sum(u0)) ## starting point of the order parameter(sum of all oscillators)
tspan = (0.0,100.0)
# K=coupling constant 
#wi=natural frequency of oscillator i 
p = (K=0.9,w2=7,w3=10,w4=13)
prob = ODEProblem(parametrized_strogatz,u0,tspan,p)
sol = solve(prob)
p1=plot(sol,vars=(0,[1,2,3,4]),labels=["Order parameter" "oscillator 1" "oscillator 2" "oscillator 3"])

anotherexample

2

Now just do a loop.

3

Hey, I am so grateful you developed this ecosystem. Many thanks.

I am afraid I do not get what you mean. Could you please explain it a little? What do I loop?

The idea is that all 800 oscillators should be instantiated simultaneously so they interact with each other from t=0

4

For anyone interested, this is the solution(I came up with) using loops as @ChrisRackauckas suggested. Thank you.

function parametrized_strogatz(du,u,p,t)
    #Coupling constant
    K=p["K"]
    ##centroid
    du[1]=0
    for i in 1:p["N"]
        du[1]+=du[i+1]
    end
    du[1]=du[1]/p["N"]
    #oscillators
    for i in 1:p["N"]
        du[i+1]=(1.0-abs2(u[i+1])+p["w$(i+1)"]*im)u[i+1]+K*(u[1]-u[i+1])
    end
end
2 Likes
5

I am not sure why you dont use a vector for the weights. Also, a named tuple for p seems cleaner.

par = (K = 1., N = 10, w = rand(10)) # sent to the following function

function parametrized_strogatz(du,u,p,t)
    #Coupling constant
    K = p.K
    N = p.N
    w = p.w

    ##centroid
    du[1]=0
    for i in 1:N
        du[1]+=du[i+1]
    end
    du[1]=du[1]/N
    #oscillators
    for i in 1:N
        du[i+1]=(1.0-abs2(u[i+1])+w[i+1]*im)*u[i+1]+K*(u[1]-u[i+1])
    end
end
2 Likes
6

I was doing more things with other functions that I did not post here, that is why. But yes, it is much cleaner the way you write it. Thank you