Abstract

This paper mainly attempts to recreate Kuramoto's model of coupled oscillators as defined in [strogatz]. This paper consists of four sections-

• §1. Describe and implement the model as in the paper

• §2. Check that the model evolves as we expect it to and make pretty animations

• §3. Show that simulated results agree with certain analytical properties

• §4. Study the robustness of the model to different distributions

How to read the paper

This paper is made using the Pluto[pluto] Package which provides an interface which can show both code and theory side by side, which will hopefully enhance the understanding of the work and also provide a reproducable implementation without having to go into the details of it. Those who wish to read the paper for the theory and results can easily skip over the code and still understand everything. Those who are attempting to recreate the work can use the code as a guide. However, in the interest of keeping the code aligned to the topic at hand, certain auxillary functions have been hidden. The content of these hidden cells can be found in the executable version of this notebook here

§0 Packages used

The Plots.jl[plots] package is used for plotting and the Statistics[stats] Package for some statistical analysis. Distributions.jl[distributions] is used to sample the . PlutoUI.jl[plutoui] is used for some UI elements such as the ToC. JLD[jld] is used for storing data between runs. Images.jl[images] is used for Gaussian smoothing. LsqFit.jl[lsqfit] is used to fit the ${\kappa }_c$ vs $\sigma$ function.

hidden code cell

§1 - The Kuramoto Model

This is a model which attempts to simulate coupled oscillators. Consider a set of $n$ oscillators $A_{i|1\le i\le n}$. Each $A_i$ has an angular position and internal angular velocity ${\theta }_i$ and ${\omega }_i$. Then, the most general coupled oscillator equation would be

$${\dot{\theta }}_j={\omega }_j+f({\theta }_i\dots {)}_{1\le i\le n}$$

If the coupling is weak, and symmetric, this can be reduced to

$${\dot{\theta }}_j={\omega }_j+\kappa \sum _{i=1}^n\sin ({\theta }_j-{\theta }_i)$$

Which can be recast into the following equation

$${\dot{\theta }}_j={\omega }_j+\kappa r\sin (\mathrm{\Psi }-{\theta }_i)$$

$$\text{where}$$

$$r{e}^{\iota \mathrm{\Psi }}=\text{mean}(e^{\iota {\theta }_i}\dots {)}_{1\le i\le n}$$

Agent

Let us define a structure called Agent which holds two variables R and V, which correspond to $\theta$ and $\omega$ respectively.

Now let's define how an agent updates, given κ, r, Ψ and dt, as defined in the equations before

Helper Functions

Evolution and Simulation

As always, we define a single step evolution along with a multistep simulate function.

The evolve function takes a list of agents and updates them according to the dynamics explained above. The simulate function sets the parameters and records certain variables for plotting and further analysis

§2 Visualizing the model

Let us first make sure that the model is running correctly. We can either visualize the positions of the oscillators with time, or we can plot $r$ vs time and make sure that the results match what is given in the paper. We will do both.

$r_{\mathrm{\infty }}$ depends on κ

Animated visualization of the dynamics

For this, we define a simulate2 function which is of the same form as the previous function, but it also records the agents for each time step. Then we plot the locations of the oscillators and the mean vector for each time step and make an animation. Due to computational limitations, only 100 particles for 10000 timesteps are taken, and animation is made for every 20th frame. Also, for reproducability, the random seed is fixed. The seed was selected so that the results were clearer to show. However all observations are observed for any number of unseeded runs. To test, simply comment out the seed line.

hidden code cell

hidden code cell

§3 Simulations agree with analysis

In section 3.3 of [strogatz], the author mentions the following observations of the simulations.

r∞ vs κ

• For all $K$ less than a certain threshold $K_c$, the oscillators act as if they were uncoupled: the phases become uniformly distributed around the circle, starting from any initial condition. Then $r(t)$ decays to a tiny jitter of size $\mathcal{O}(N^{-1/2})$

• When $K>K_c$, this incoherent state becomes unstable and $r(t)$ grows exponentially, reflecting the nucleation of a small cluster of oscillators that are mutually synchronized, thereby generating a collective oscillation. Eventually $r(t)$ saturates at some level $r_{\mathrm{\infty }}<1$, though still with $\mathcal{O}(N^{-1/2})$ fluctuations.

Both of these claims are shown to be true in the first plot. For a "high" $K(=7)$, $r$ settles with jitter at 1, and for a "low" $K(=1)$, $r$ settles with jitter at 0.

Individual Oscillators

• The population splits into two groups: the oscillators near the center of the frequency distribution lock together at the mean frequency and co-rotate with the average phase ψ(t), those in the tails run near their natural frequencies and drift relative to the synchronized cluster.

• This mixed state is often called partially synchronized. With further increases in K, more and more oscillators are recruited into the synchronized cluster, and $r_{\mathrm{\infty }}$ tends to 1

Both of these claims are also shown to be true in the animation and in the plots of the effective velocities. In the animation, we clearly see that most of the oscillators are in a clump, but there are a few that drift away from the mean and stay out of the cluster. In the velcity plot, this is almost clearly evident. Also, if looked at carefully, it shows that the particles which are not in sync are the ones which where far away from the mean in the initial distribution

κc vs g(ω)

The paper also analyses ${\kappa }_c$ as a function of $g(\omega )$. While Kuramoto and Nishikawa had previously assumed that the drifters would not have any effect on the coherence of the oscillators, they later recanted their stance and assumed that the drifters did in fact have an effect on the evolution of the coherence. Probably if they had the luxury of fast computers and easy computer languages, they would have seen that the drifters did in fact have an effect on $r(t)$. In the animation above, there are small but fast changes in the angle and the amplitude of the mean vector, and these are directly related to the movement of the drifters. Obviously, a simulation of only 100 oscillators is obviously not in the thermodynamic limit and the effect on $r$ is more apparent.

None the less, we can be reasonably sure that Kuramoto and Nishikawa's second attempt (§6.2 in [strogatz]) was more accurate, and by extension so was Strogatz's own analysis. From their analysis,

$${\kappa }_c=\frac{2}{\pi g(0)}$$

In the case of the gaussian distribution,

$$g(0)=\frac{1}{\sigma \sqrt{2\pi }}$$

Hence,

$${\kappa }_c=\frac{2\sigma \sqrt{2\pi }}{\pi }=\sqrt{\frac{8}{\pi }}\sigma$$

Computational implementation

We define a function simulate_for_sigma which searches by bisection for $\kappa$s which are useful for our plot. We limit the difference of $\kappa$ and difference in $r_{\mathrm{\infty }}$ for the sake of computational time. $r_{\mathrm{\infty }}$ is evaluated by taking a mean over the last $10\mathrm{\%}$ timesteps.

In the ideal case, we'd average $r_{\mathrm{\infty }}$ over many initial conditions to make sure that there are infact no artifacts of the initial conditions. This is especially true, near the bifurcation, where the transition times can be very long. However, we take the easier path by simply box-smoothing over the data (a tool often employed by statisticians when they get noisy data).

Run new simulation?

Stability of κc over σ

Using the powerful computers at our disposal, we can run the simulation on multiple $\kappa$s and $\sigma$s and find the $r_{\mathrm{\infty }}$ for each of these. Then, we can analyse the results to find ${\kappa }_c$ for each $\sigma$ and fit it to a straight line and see if the slope is equal to $\sqrt{\frac{8}{\pi }}$. Like before, we smooth our data, but this time we convolve the data with a 2dimensional Gaussian blur filter.

We measure the points of phase transition by taking a derivative along $\kappa$, and picking points where the derivative is above a certain threshold(=0.09).

Code to generate data

function graph(κs::Array{Float64}, σs::Array{Float64})
    arr = zeros(Float64, length(κs), length(σs))
    for (i, κ) in enumerate(κs), (j, σ) in enumerate(σs)
        arr[i,j] = mean(
                abs, 
                @view simulate(Int(1e5), 1000, σ, κ, 1e-4)[end - 1000:end]
        )
    end
    arr
end
κs, σs = collect(0.0:0.1:2.5), collect(0.0:0.1:2.5)
zs = graph(κs, σs)
save("save.jld", "zs", zs, "ks", κs, "ss", σs)

Interactive plot below

Interactive plot below

$${\kappa }_c=1.56\cdot \sigma$$

Observations

The single run and the extended analysis both show the same result that the ${\kappa }_c$ follows the analysis in the paper closely, however, there is a small disparity. Any finite simulation will inevitably not have the sharp nature of the bifurcation. In §4, we run the same simulations for other distributions, and see the effect of the kind of distribution on this.

§4 Robustness of results to distributions

As mentioned in the previous subsection, let's run the same simulation with some other distributions

Gaussian Distribution

We've already run this before, but let's do it again.

$$g(\omega )=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp [-\frac{{\omega }^2}{2{\sigma }^2}];g(0)=\frac{1}{\sqrt{2\pi {\sigma }^2}};{\kappa }_c(3)=\sqrt{\frac{8}{\pi }}\ast 3=1.596$$

Lorenz-Cauchy distributions

$$g(\omega )=\frac{1}{\pi \sigma (1+{\left(\frac{\omega }{\sigma }\right)}^2)};g(0)=\frac{1}{\pi \sigma };{\kappa }_c(3)=6$$

Reflected exponential decay

$$g(\omega )=\frac{1}{2a}\exp (-\frac{|x|}{a});g(0)=\frac{1}{2a};{\kappa }_c(3)=\frac{4\pi }{3}=4.189$$

Observations

Like the Gaussian, the other distributions also have a slight discrepancy from the expected value. However by a visual analysis, we can say that the theoretical results are close to what we would get in the infinite limit. For a more accurate analysis, one could fit the function as mentioned in the paper to the data, and then look at the value of ${\kappa }_c$.

Another observation that we can see is that distributions that are tail heavy are much more prone to variations in the partially coherent part of the graph. In order to reduce this, we average results over 3 runs with random initial conditions for the Lorenz Cauchy distribution

References