Simulates the N-oscillator Kuramoto model

kuramoto.py

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

N = 5
K = 12
ω = np.random.normal(0, 1, N)
θ_0 = np.random.uniform(0, 2*np.pi, N)

def f(t, θ):
  return [ω[i] + K/N * sum(np.sin(θ[i] - θ[j]) for j in range(N)) for i in range(N)]

t_max = 10
dt = 0.1
sol = integrate.solve_ivp(f, t_span=[0,t_max], y0=θ_0, t_eval=np.arange(0,t_max,dt))

for i in range(N):
  plt.plot(sol.t, sol.y[i] % 2*np.pi)
plt.xlabel("t")
plt.ylabel("θ")
plt.show()