The Kuramoto synchronization model for coupled oscillators¶
Our goal is to explore a model that is used to described synchronization of oscillators.
The Kuramoto model is a set of $N$ interacting oscillators. Each oscillator is described with an angle $\theta_i$. The indices range from 0, .... $N-1$.
Each oscillator has an intrinsic frequency or phase velocity $\omega_i$.
In the absence of interactions between oscillators $$ \frac{d\theta_i}{dt} = \omega_i$$
The phases continually advance, as you would expect for the phase of an isolated harmonic oscillator.
With interactions $$\frac{d \theta_i}{dt} = \omega_i + \sum_{j=0}^{N-1} K_{ij} \sin (\theta_j- \theta_i) $$
If the interaction strengths $K_{ij}>0$ are high enough the oscillators will synchronize.
Notice that all the oscillators interact with every other oscillator. The model is more appropriate for a bunch of fireflies than a chain of oscillators with nearest neighbor interactions.
A simple version of the Kuramoto coupled oscillator model¶
For our simple exploration we set $\omega_i=1$ for all oscillators and we set $K_{ij} = K>0$ the same for each pair of oscillators.
This gives $$\frac{d \theta_i}{dt} = 1 + \sum_{j=0}^{N-1} K \sin (\theta_j- \theta_i) $$
We start with an array of phases $\theta_i$. We use the index $n$ to refer to the time step of our integrator. From the phases $\theta_i^n$, at time step $n$ we want to compute the phases at the next time step. If the timestep is $\Delta t$
$$\frac{d \theta^n_i}{dt} \sim \frac{ \theta^{n+1}_i - \theta^n_i}{\Delta t} $$This gives $$ \theta^{n+1}_i = \theta_i^n + \Delta t \frac{d \theta^n_i}{dt} $$ With the equations of motion this gives
$$ \theta^{n+1}_i = \theta_i^n + \Delta t \left(1 + \sum_j K \sin (\theta_j^n- \theta_i^n)\right) $$We can use this to compute the new array of phases from the previous one.
Our goal¶
Our goal is to write a routine that takes
as arguments, an array of $N$
phases each between 0 and $2\pi$,
the interaction strength $K$ and a time step interval $\Delta t$.
The routine should return the new phase array after a time $\Delta t$!
It would also be useful to have a routine that returns arrays at different timesteps or numbers of iterations of the previous routine.
Some notes:
You need to keep the phases within $[0,2\pi)$. Take the modulo $2 \pi$ after you compute the new phases.
Keep the timestep $\Delta t \ll 1$.
initial conditions¶
I would initialize the array with randomly chosen phases. For example with $N$ the number of oscillators
theta_arr_init = random.random(N)*2*np.pigives initial phases randomly chosen in $[0,2\pi)$
I found that it can be faster to compute the new array after a single timestep using the function numpy.roll()
https://numpy.org/doc/stable/reference/generated/numpy.roll.html
For example
dtheta_arr = np.zeros(N) + 1 # intrinsic phase velocity is 1 for all oscillators
for j in range(N): # loop over all possible shifts
thetaj = np.roll(theta_arr,j) # shift the array by j
dtheta_arr += K*np.sin(thetaj - theta_arr)
new_theta_arr = theta_arr + dtheta_arr*dtwhat to plot?¶
We can fill a 2d array with phase array outputs. One dimension of the 2d array is $N$, the number of oscillators. The other dimension of the 2d array is the number of time steps.
We could then plot this 2D array as an image. We should be able to tell if all the oscillators synchronize and start to have the same phases as each other.
As all oscillators interact with all the other oscillators there is no prefered order. We could plot all the phases as dots on the unit circle and animate the plot.
You can plot an order parameter as a function of time.
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
#Here is an example of using roll to shift arrays
xarr = np.arange(10)
plt.plot(xarr,'ro',label='0')
xarr_roll_1 = np.roll(xarr,1)
plt.plot(xarr_roll_1,'go',label='1')
xarr_roll_2 = np.roll(xarr,2)
plt.plot(xarr_roll_2,'bo',label='1')
xarr_roll_m1 = np.roll(xarr,-1)
plt.plot(xarr_roll_m1,'co',label='-1')
plt.legend()
