<jupyter_start><jupyter_text>KuramotoWe already found that complex systems show signs of feedback dynamics. In syncronisation one explores how the nonlinear feedback can be used to syncronise a population of weakly coupled entities. In this notebook we are going to explore the Kuramoto model which shows syncronisation behaviour but still exactly solvable under specific restrictions.<jupyter_code>from plotly import offline as py
from plotly import graph_objs as go

py.init_notebook_mode(connected=True)<jupyter_output><empty_output><jupyter_text>The Kuramoto model describes a population of $N$ weakly couped oscillators. Each oscillator evolves according to, $$\dot{\theta}(t)=\omega_i+\frac{K}{N}\sum^N_{j=1}\sin\left(\theta_j-\theta_i\right),$$ wherein $\omega_i$ is the fundamental frequency of the oscillator and $K$ denotes the mean field coupling between the oscillators. We take the fundamental frequency to follow a standard Cauchy (Lorentzian) probability distribution $\mathcal{C}(x_0=0,\gamma=1)$. The probability density function is provided by$$f(\omega)=\frac{1}{\pi}\frac{\gamma}{\omega^2+1}.$$<jupyter_code>from scipy.stats import cauchy

def kuramoto(t, K, N, dh=1e-1):
    omega = cauchy.rvs(size=N)
    times = np.arange(0, t, dh)
    theta = np.zeros((len(times), N))
    
    for i in range(1, len(times)):
        for j in range(N):
            dtheta = omega[j] + K * np.mean(np.sin(theta[i-1] - theta[i-1][j]))
            
            theta[i][j] = theta[i-1][j] + dtheta * dh
            
    return times, theta, np.exp(1j*theta).mean(axis=-1), omega<jupyter_output><empty_output><jupyter_text>We integrate the oscillator using the Euler method.<jupyter_code>time, theta, order, _ = kuramoto(100, 0, 100)

r = np.abs(order)
phi = np.degrees(np.arctan2(order.imag, order.real))

py.iplot([
    go.Scatterpolar(r=r, theta=phi, mode='markers', marker=dict(color=time, colorbar=dict(title='time'), colorscale='Viridis')),
])<jupyter_output><empty_output><jupyter_text>We note that the oscillator start in a syncronized state from which they move into a state where the oscillator are out of sync and their global phase moves moves around.<jupyter_code>time, theta, order, _ = kuramoto(100, 5, 100)

r = np.abs(order)
phi = np.degrees(np.arctan2(order.imag, order.real))

py.iplot([
    go.Scatterpolar(r=r, theta=phi, mode='markers', marker=dict(color=time, colorbar=dict(title='time'), colorscale='Viridis')),
])<jupyter_output><empty_output><jupyter_text>For intermediate coupling we find the oscilator to transition into a state where they are in a semi ordered state.<jupyter_code>time, theta, order, _ = kuramoto(100, 10, 100)

r = np.abs(order)
phi = np.degrees(np.arctan2(order.imag, order.real))

py.iplot([
    go.Scatterpolar(r=r, theta=phi, mode='markers', marker=dict(color=time, colorbar=dict(title='time'), colorscale='Viridis')),
])<jupyter_output><empty_output><jupyter_text>For larger coupling $K/N=0.1$ the oscillator approach an ordered state.<jupyter_code>time1, theta1, order1, _ = kuramoto(100, 5, 10)
time2, theta2, order2, _ = kuramoto(100, 5, 100)
time3, theta3, order3, _ = kuramoto(100, 5, 1000)

r1 = np.abs(order1)
r2 = np.abs(order2)
r3 = np.abs(order3)

phi1 = np.degrees(np.arctan2(order1.imag, order1.real))
phi2 = np.degrees(np.arctan2(order2.imag, order2.real))
phi3 = np.degrees(np.arctan2(order3.imag, order3.real))

py.iplot([
    go.Scatterpolar(r=r1, theta=phi1, mode='markers', name='N = 10'),
    go.Scatterpolar(r=r2, theta=phi2, mode='markers', name='N = 100'),
    go.Scatterpolar(r=r3, theta=phi3, mode='markers', name='N = 1000'),
])<jupyter_output><empty_output><jupyter_text>We note that the spread of the order parameter is decreases with the number of oscillators. Furthermore it appears that the global phase drifts more slowly for large number of oscillators.<jupyter_code>K = np.linspace(0, 10)
r = []

for (i, k) in enumerate(K):
    time, theta, order, _ = kuramoto(100, k, 100)
    
    r.append(np.abs(order[-1]))
layout = go.Layout(
    title='Order parameter',
    xaxis=dict(title='K'),
    yaxis=dict(title='r(K)'),
)

py.iplot([
    go.Scatter(x=K, y=r, mode='lines+markers', name='experiment'),
    go.Scatter(x=K, y=np.sqrt(1-2/K), mode='lines+markers', name='prediction'),
])<jupyter_output><empty_output><jupyter_text>For Cauchy distributed intrinsic oscillator frequencies $\omega_i$ the theory predicts that the order parameter converges to, $$r(K)=\sqrt{1-\frac{K^*}{K}},$$wherein $K^*$ denotes the critical coupling value which is linked to the probability density function of the intrinsic frequencies $f(\omega)$ through,$$K^*=\frac{2}{\pi f(0)}=2,$$with the last equal being true for our Lorentzian distribution.<jupyter_code>time, theta, _, omega = kuramoto(100, 5, 100)

order = np.exp(1j * theta)

r = np.abs(order)
phi = np.degrees(np.arctan2(order.imag, order.real))

py.iplot([
    go.Scatterpolar(r=r[-1], theta=phi[-1], mode='markers', marker=dict(color=omega, colorbar=dict(title='ω'), colorscale='Viridis'))
])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1462787015043547, 'fraction_numerical': 0.02632620744259699, 'mean_word_length': 5.436942675159235, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '193148', 'license': 'Apache-2.0', 'n_tokens_mistral': 1891, 'n_tokens_neox': 1736, 'n_words': 456, 'path': 'notebooks/syncronisation/Kuramoto.ipynb', 'repo_name': 'bodokaiser/complex-systems', 'size': '5052'}