The Kuramoto synchronization model

# import numpy and a random number generator
import numpy as np
from numpy import random
# make it possible to plot figures within the notebook
%matplotlib inline
from matplotlib import pyplot as plt
# matplotlib allows us to make figures

twopi = 2*np.pi  # useful global variable

# update a time step!
# arguments: 
#     a theta_arr array of phases for N oscillators
#     K interaction strength K
#     omega_arr can be an array or a single number, lists intrinsic frequencies
#      timestep dt
# returns: a new theta array after integrating 1 timestep
def new_theta(theta_arr,omega_arr,K,dt):
    N = len(theta_arr)
    new_theta_arr = np.copy(theta_arr)
    dtheta = np.zeros(len(theta_arr))
    # looping over all pairs of oscillators
    for j in range(1,N):  # for each possible shift
        shift_theta_arr = np.roll(theta_arr,j)  # periodic boundary
        dtheta += K*np.sin(shift_theta_arr - theta_arr)
        # adding to whole array at once
        # this is faster than doing a double loop
    dtheta += omega_arr  # these are the intrinsic frequencies  
    # instrisic frequencies could be an array of different frequencies
    new_theta_arr += dt*dtheta
    return new_theta_arr%twopi  # modulo 2 pi

# carry out nsteps timesteps of the Kuramoto model
# returns: new array of phases for each oscillator
# arguments: 
#    theta_arr has phases in it
#    K is the interaction strength
#    dt is the time step
#    nsteps  how many timesteps to do
#    omega_arr an array of intrinsic frequencies
def do_nsteps(theta_arr,omega_arr,K,dt,nsteps):
    new_theta_arr = np.copy(theta_arr)
    for i in range(nsteps):   # iterate nsteps
        new_theta_arr = new_theta(new_theta_arr,omega_arr,K,dt)
    return new_theta_arr  # return new array!
        
# fill a 2d array with timesteps
# first index gives which oscillator
# second index gives which timestep
# returns: 
#    the 2d array, zarr
# arguments: 
#    theta_arr has initial conditions (phases)
#    K is the interaction strength
#    dt is the time step
#    nsteps, how many timesteps to do
#    omega_arr an array of intrinsic frequencies
def fill2d(theta_arr,omega_arr,K,dt,nsteps):
    zarr = np.zeros((N,nsteps))   # our 2d array
    zarr[:,0] = theta_arr  # initial conditions
    new_theta_arr = theta_arr
    for i in range(1,nsteps):
        new_theta_arr = new_theta(new_theta_arr,omega_arr,K,dt)
        zarr[:,i] = new_theta_arr
     
    return zarr
    

# just plotting phases as a function of index $j$
N=20
theta_arr = 2*np.pi*np.random.random(N) # random initial  conditions
K=0.1; dt=0.1
nsteps=15
omega_arr = np.ones(N)  # intrinsic frequency array

new_theta_arr=do_nsteps(theta_arr,omega_arr,K,dt,nsteps)
plt.plot(new_theta_arr,'ro')
plt.ylim((0,twopi))
new_theta_arr=do_nsteps(new_theta_arr,omega_arr,K,dt,nsteps)
plt.plot(new_theta_arr,'go')
new_theta_arr=do_nsteps(new_theta_arr,omega_arr,K,dt,nsteps)
plt.plot(new_theta_arr,'bo')
plt.xlabel('index')
plt.ylabel('phase')

Text(0, 0.5, 'phase')

# make a 2d array with timesteps
N=20
theta_arr = 2*np.pi*np.random.random(N)
K=0.01; dt=0.1
nsteps=200
omega_arr = np.ones(N) 
zarr = fill2d(theta_arr,omega_arr,\
              K,dt,nsteps)  # fill a 2d array [i,j] #i is oscillator, j is timestep
tmax = nsteps*dt
extent = ((0,tmax,0,N))
fig,ax = plt.subplots(1,1,figsize=(5,3))
ax.imshow(zarr,extent=extent)
ax.set_xlabel('time')
ax.set_ylabel('oscillator')
#plt.colorbar(fraction=0.02)

Text(0, 0.5, 'oscillator')

from matplotlib import animation, rc
from IPython.display import HTML

# First set up the figure, the axis, and the plot element we want to animate
fig, ax = plt.subplots()

rmax=1.1
ax.set_xlim((-rmax, rmax))
ax.set_ylim((-rmax, rmax))
ax.set_aspect('equal')

# this is where plotting object is defined!
dots, = ax.plot([],[],'ro', lw=2,ms=6,alpha=0.5)

# initialization function: plot the background of each frame
# line is a line2D object that we have previously defined globally 
def init():
    dots.set_data([], [])
    animate_fun(0)
    return (dots,)   # returns a tuple with the object *line* in it

# animation function. This is called sequentially with integer i increasing
# this routine expects that a global object line exists already
def animate_fun(i):
    theta_arr = np.squeeze(zarr[:,i])
    x = np.cos(theta_arr)
    y = np.sin(theta_arr)
    dots.set_data(x, y)
    return (dots,)

# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate_fun, init_func=init,
                               frames=200, interval=20, blit=True)

HTML(anim.to_jshtml()) 

Once Loop Reflect