Kuramoto Assignment (Graded)

Module 2: Synchronization and Collective Behaviour

This is a graded assignment, released in Session 12. It is worth 10% of your final grade.

Use of generative AI

Following the syllabus AI policy, you may not use generative AI tools to solve this assignment. This covers the code, the analysis and the written report.

Your goal is to review the theory of synchronization of weakly coupled oscillators, study the Kuramoto model numerically, and then simulate the Millennium Bridge crowd-bridge synchrony model (Kuramoto 2003; Strogatz 2000; Pikovsky, Rosenblum, and Kurths 2001; Strogatz et al. 2005; Eckhardt et al. 2007). The lab pages of Session 12 give you the building blocks. The theory is in Critical Coupling and The Millennium Bridge.

Kuramoto Model

Start from the basic Kuramoto model with all-to-all symmetric coupling:

\[ \dot\theta_i = \omega_i + \frac{K}{N}\sum_{j=1}^N \sin(\theta_j-\theta_i), \qquad i=1,\dots,N. \tag{1}\]

Here $\theta_i \in [0,2\pi)$ is the phase of oscillator $i$, $K>0$ is the coupling constant, and the natural frequencies $\omega_i$ are sampled from a unimodal symmetric distribution $g(\omega)$.

Numerical Integration

Write a Python script that integrates the Kuramoto equations numerically. Start with the following values:

$N=100$ oscillators.
$g(\omega)$ equal to the standard normal distribution.
$K=7$.

For the initial condition, write a function that assigns phases either randomly over the whole circle, the disperse case, or over a small arc, the concentrated case. For the numerical integration, you may write your own scheme, for instance explicit Euler, or use scipy.integrate. In either case, produce the following outputs:

A graph of $r(t)$, the order parameter, versus time. In the subcritical case $K<K_c$, the asymptotic value $r_\infty$ should be approximately zero up to fluctuations, while for $K>K_c$ you should observe $r_\infty \in (0,1)$ for large $t$. The amplitude of the fluctuations should decrease as $N$ grows, typically like $\mathcal{O}(N^{-1/2})$.
An animation of the oscillators moving on the unit circle. You may also plot the centroid $re^{i\psi}$.

Critical Coupling

Review Critical Coupling and Sections 3 and 4 of Strogatz (2000) until you are comfortable with the derivation of the critical coupling $K_c$. For more detail, see Pikovsky, Rosenblum, and Kurths (2001). Then compare the theoretical prediction with the value obtained numerically.

Sample 20 values of $K$ between $K_{\min}$ and $K_{\max}$ such that $K_{\min}<K_c<K_{\max}$. A reasonable choice is $K_{\min}=K_c/3$ and $K_{\max}=3K_c$.
For each sampled value of $K$, estimate $r_\infty$ by making sure that $r(t)$ has reached its saturation regime, then average over a tail window to reduce fluctuations. A reasonable choice is to average over the last 300 time steps. Keep both the mean and standard deviation.
Plot $r_\infty(K)$ with error bars and identify the approximate value of $K_c$, comparing it with the theoretical prediction.
Validate your results with different choices of $g(\omega)$, for instance normal distributions with different variance, and check whether the theoretical prediction for $K_c$ remains consistent.
Plot the smooth theoretical curve defined implicitly by the self-consistency equation (4.5) of Strogatz (2000) together with your numerical points:

\[ 1 = K\int_{-\pi/2}^{\pi/2} \cos^2\theta\, g\big(Kr\sin\theta\big)\,d\theta. \tag{2}\]

If you are not able to solve the implicit equation numerically for a general distribution $g(\omega)$, include at least the Cauchy case, where the formula is explicit:

\[ g(\omega) = \frac{1}{\pi(1+\omega^2)}, \qquad r_\infty(K) = \sqrt{1-\frac{K_c}{K}}. \tag{3}\]

Check that your numerical points approach the theoretical curve more closely as $N$ becomes large.

Millennium Bridge Crowd Synchrony

Read The Millennium Bridge first. The model couples one lateral bridge mode to $N$ pedestrian phase oscillators (Strogatz et al. 2005):

\[ M\ddot X + B\dot X + KX = G\sum_{i=1}^N \sin\Theta_i, \qquad \dot\Theta_i = \Omega_i + C\,A\sin\big(\Psi-\Theta_i+\alpha\big), \tag{4}\]

with $\Omega_0 = \sqrt{K/M}$ and the bridge amplitude and phase given by

\[ A = \sqrt{X^2 + \big(\dot X/\Omega_0\big)^2}, \qquad \Psi = \operatorname{atan2}\big(X,\ \dot X/\Omega_0\big). \tag{5}\]

Here $K$ is the bridge stiffness, not the Kuramoto coupling. Crowd synchrony is measured with $R e^{i\Phi} = N^{-1}\sum_j e^{i\Theta_j}$.

For $\alpha = \pi/2$ and a frequency density $P(\Omega)$ symmetric about $\Omega_0$, the critical crowd size is

\[ N_c = \frac{2B\,\Omega_0}{\pi\,G\,C\,P(\Omega_0)}. \tag{6}\]

If $P$ is normal with mean $\Omega_0$ and standard deviation $\sigma$, then $P(\Omega_0) = 1/(\sqrt{2\pi}\,\sigma)$.

Required Numerical Experiment

Implement a numerical solver for the coupled bridge-pedestrian system. You may use scipy.integrate.solve_ivp, which is recommended here, or another standard ODE solver such as RK4. Use random initial phases $\Theta_i(0) \sim \mathrm{Unif}[0,2\pi)$ and small initial bridge motion, for instance $X(0) \approx 0$ and $\dot X(0) \approx 0$.

Perform the following controlled crowd-ramp experiment:

Start with $N=N_0$ pedestrians. Every $\Delta T$ seconds, increase the crowd size by $\Delta N$, adding new pedestrians with fresh initial phases and natural frequencies, until reaching $N_{\max}$.
Compute the theoretical value of $N_c$ using your chosen parameters and frequency distribution.
Define the first time at which the crowd reaches the theoretical threshold,

\[ t_c = \inf\{t : N(t) \ge N_c\}. \tag{7}\]

Produce three stacked plots, and in each one draw a vertical dashed line at $t=t_c$ labeled “$N=N_c$”:
- $N(t)$ versus $t$, the crowd-size staircase.
- $A(t)$ versus $t$, the bridge wobble amplitude.
- $R(t)$ versus $t$, the degree of synchronization.
Briefly describe what happens before and after $t_c$, and comment on whether the observed onset of wobbling and synchrony is close to the theoretical prediction.

Reference code is available in amlab/odes_coupled/bridge.py. It uses a different, nondimensional formulation of the bridge-pedestrian coupling, so its parameters do not map one to one onto the equations above. Use it for ideas, but do not copy it.

What to submit

Your code, which must run from start to finish and reproduce every figure (fix your random seeds).
A short report with the figures requested above and your interpretation.

References

Eckhardt, Bruno, Edward Ott, Steven H. Strogatz, Daniel M. Abrams, and Alan McRobie. 2007. “Modeling Walker Synchronization on the Millennium Bridge.” Physical Review E 75: 021110.

Kuramoto, Yoshiki. 2003. Chemical Oscillations, Waves, and Turbulence. Dover Publications.

Pikovsky, Arkady, Michael Rosenblum, and J"urgen Kurths. 2001. Synchronization: A Universal Concept in Nonlinear Sciences. Cambridge University Press.

Strogatz, Steven H. 2000. “From Kuramoto to Crawford: Exploring the Onset of Synchronization in Populations of Coupled Oscillators.” Physica D: Nonlinear Phenomena 143: 1–20.

Strogatz, Steven H., Daniel M. Abrams, Alan McRobie, Bruno Eckhardt, and Edward Ott. 2005. “Crowd Synchrony on the Millennium Bridge.” Nature 438: 43–44.

--- title: "Kuramoto Assignment (Graded)" subtitle: "Module 2: Synchronization and Collective Behaviour" format: html --- This is a **graded assignment**, released in Session 12. It is worth 10% of your final grade.  ::: {.callout-warning} ## Use of generative AI Following the syllabus AI policy, you may not use generative AI tools to solve this assignment. This covers the code, the analysis and the written report. ::: Your goal is to review the theory of synchronization of weakly coupled oscillators, study the Kuramoto model numerically, and then simulate the Millennium Bridge crowd-bridge synchrony model [@kuramoto2003chemical; @strogatz2000kuramoto; @pikovsky2001synchronization; @strogatz2005bridge; @eckhardt2007modeling]. The lab pages of Session 12 give you the building blocks. The theory is in [Critical Coupling](critical-coupling.qmd) and [The Millennium Bridge](millennium-bridge.qmd). ## Kuramoto Model Start from the basic Kuramoto model with all-to-all symmetric coupling: $$ \dot\theta_i = \omega_i + \frac{K}{N}\sum_{j=1}^N \sin(\theta_j-\theta_i), \qquad i=1,\dots,N. $$ {#eq-ode-coupled-assignment-1} Here $\theta_i \in [0,2\pi)$ is the phase of oscillator $i$, $K>0$ is the coupling constant, and the natural frequencies $\omega_i$ are sampled from a unimodal symmetric distribution $g(\omega)$. ## Numerical Integration Write a Python script that integrates the Kuramoto equations numerically. Start with the following values: - $N=100$ oscillators. - $g(\omega)$ equal to the standard normal distribution. - $K=7$. For the initial condition, write a function that assigns phases either randomly over the whole circle, the disperse case, or over a small arc, the concentrated case. For the numerical integration, you may write your own scheme, for instance explicit Euler, or use `scipy.integrate`. In either case, produce the following outputs: 1. A graph of $r(t)$, the order parameter, versus time. In the subcritical case $K<K_c$, the asymptotic value $r_\infty$ should be approximately zero up to fluctuations, while for $K>K_c$ you should observe $r_\infty \in (0,1)$ for large $t$. The amplitude of the fluctuations should decrease as $N$ grows, typically like $\mathcal{O}(N^{-1/2})$. 2. An animation of the oscillators moving on the unit circle. You may also plot the centroid $re^{i\psi}$. ## Critical Coupling Review [Critical Coupling](critical-coupling.qmd) and Sections 3 and 4 of @strogatz2000kuramoto until you are comfortable with the derivation of the critical coupling $K_c$. For more detail, see @pikovsky2001synchronization. Then compare the theoretical prediction with the value obtained numerically. 1. Sample 20 values of $K$ between $K_{\min}$ and $K_{\max}$ such that $K_{\min}<K_c<K_{\max}$. A reasonable choice is $K_{\min}=K_c/3$ and $K_{\max}=3K_c$. 2. For each sampled value of $K$, estimate $r_\infty$ by making sure that $r(t)$ has reached its saturation regime, then average over a tail window to reduce fluctuations. A reasonable choice is to average over the last 300 time steps. Keep both the mean and standard deviation. 3. Plot $r_\infty(K)$ with error bars and identify the approximate value of $K_c$, comparing it with the theoretical prediction. 4. Validate your results with different choices of $g(\omega)$, for instance normal distributions with different variance, and check whether the theoretical prediction for $K_c$ remains consistent. 5. Plot the smooth theoretical curve defined implicitly by the self-consistency equation (4.5) of @strogatz2000kuramoto together with your numerical points: $$ 1 = K\int_{-\pi/2}^{\pi/2} \cos^2\theta\, g\big(Kr\sin\theta\big)\,d\theta. $$ {#eq-ode-coupled-assignment-2} 6. If you are not able to solve the implicit equation numerically for a general distribution $g(\omega)$, include at least the Cauchy case, where the formula is explicit: $$ g(\omega) = \frac{1}{\pi(1+\omega^2)}, \qquad r_\infty(K) = \sqrt{1-\frac{K_c}{K}}. $$ {#eq-ode-coupled-assignment-3} Check that your numerical points approach the theoretical curve more closely as $N$ becomes large. ## Millennium Bridge Crowd Synchrony Read [The Millennium Bridge](millennium-bridge.qmd) first. The model couples one lateral bridge mode to $N$ pedestrian phase oscillators [@strogatz2005bridge]: $$ M\ddot X + B\dot X + KX = G\sum_{i=1}^N \sin\Theta_i, \qquad \dot\Theta_i = \Omega_i + C\,A\sin\big(\Psi-\Theta_i+\alpha\big), $$ {#eq-ode-coupled-assignment-4} with $\Omega_0 = \sqrt{K/M}$ and the bridge amplitude and phase given by $$ A = \sqrt{X^2 + \big(\dot X/\Omega_0\big)^2}, \qquad \Psi = \operatorname{atan2}\big(X,\ \dot X/\Omega_0\big). $$ {#eq-ode-coupled-assignment-7} Here $K$ is the bridge stiffness, not the Kuramoto coupling. Crowd synchrony is measured with $R e^{i\Phi} = N^{-1}\sum_j e^{i\Theta_j}$. For $\alpha = \pi/2$ and a frequency density $P(\Omega)$ symmetric about $\Omega_0$, the critical crowd size is $$ N_c = \frac{2B\,\Omega_0}{\pi\,G\,C\,P(\Omega_0)}. $$ {#eq-ode-coupled-assignment-9} If $P$ is normal with mean $\Omega_0$ and standard deviation $\sigma$, then $P(\Omega_0) = 1/(\sqrt{2\pi}\,\sigma)$. ## Required Numerical Experiment Implement a numerical solver for the coupled bridge-pedestrian system. You may use `scipy.integrate.solve_ivp`, which is recommended here, or another standard ODE solver such as RK4. Use random initial phases $\Theta_i(0) \sim \mathrm{Unif}[0,2\pi)$ and small initial bridge motion, for instance $X(0) \approx 0$ and $\dot X(0) \approx 0$. Perform the following controlled crowd-ramp experiment: 1. Start with $N=N_0$ pedestrians. Every $\Delta T$ seconds, increase the crowd size by $\Delta N$, adding new pedestrians with fresh initial phases and natural frequencies, until reaching $N_{\max}$. 2. Compute the theoretical value of $N_c$ using your chosen parameters and frequency distribution. 3. Define the first time at which the crowd reaches the theoretical threshold, $$ t_c = \inf\{t : N(t) \ge N_c\}. $$ {#eq-ode-coupled-assignment-11} 4. Produce three stacked plots, and in each one draw a vertical dashed line at $t=t_c$ labeled "$N=N_c$": - $N(t)$ versus $t$, the crowd-size staircase. - $A(t)$ versus $t$, the bridge wobble amplitude. - $R(t)$ versus $t$, the degree of synchronization. 5. Briefly describe what happens before and after $t_c$, and comment on whether the observed onset of wobbling and synchrony is close to the theoretical prediction. Reference code is available in [amlab/odes_coupled/bridge.py](../../amlab/odes_coupled/bridge.py). It uses a different, nondimensional formulation of the bridge-pedestrian coupling, so its parameters do not map one to one onto the equations above. Use it for ideas, but do not copy it. ## What to submit - Your code, which must run from start to finish and reproduce every figure (fix your random seeds). - A short report with the figures requested above and your interpretation.