Kuramoto
Kuramoto oscillators coupled all to all, theta_i' = omega_i + (K / N) sum_j sin(theta_j - theta_i).
- Name:
kuramoto - Field:
physics/synchronization - Formalism: ordinary differential equations
- References: Kuramoto 1984; Datseris and Parlitz 2022; Gros 2015; Jensen 2023
Factory
@register_model(
field='physics/synchronization',
tags=('synchronization', 'oscillators', 'phase', 'mean field'),
references=('Kuramoto 1984', 'Datseris and Parlitz 2022', 'Gros 2015', 'Jensen 2023'),
defaults={'until': 50.0, 'dt': 0.5},
observables=lambda system: _order_parameter('order', 'oscillators'),
)
def kuramoto(n_oscillators: int = 50, *, K: float = 4.0, sigma: float = 1.0, seed: int = 0) -> System:
"""Return Kuramoto oscillators coupled all to all, `theta_i' = omega_i + (K / N) sum_j sin(theta_j - theta_i)`.
Every oscillator reads the group means of `cos theta` and `sin theta`, so the coupling costs time linear in the
number of oscillators. The natural frequencies are drawn from a normal distribution of width `sigma`, and the
phases uniformly from the circle. Above the critical coupling `K_c = 2 sigma sqrt(2 / pi)` a synchronized cluster
forms, and the order parameter `r`, the modulus of the mean of `exp(i theta)`, grows from about `1 / sqrt(N)`
towards 1.
Args:
n_oscillators:
The number of oscillators.
K:
The coupling strength.
sigma:
The standard deviation of the natural frequencies.
seed:
The seed of the natural frequencies and the initial phases.
Returns:
system:
The model, with the phases `oscillators.theta` and the order parameter `order`.
Examples:
>>> from skcomplex.models.physics.synchronization import kuramoto
>>> system = kuramoto(10)
>>> system.elements[0].size, system.space
(10, None)
"""
rng = np.random.default_rng(seed)
oscillators = GroupElement(
'oscillators',
element=PhaseOscillatorElement,
size=n_oscillators,
values={'theta': _spread_phases(n_oscillators, seed), 'omega': rng.normal(0.0, sigma, n_oscillators)},
)
coupling = KuramotoCouplingInteraction(
'coupling',
target='oscillators',
K=K,
mode='group',
reduce='mean',
include_self=True,
)
return System(
elements=[oscillators],
interactions=[coupling],
)
Components
class PhaseOscillatorElement(Element):
"""Phase oscillator, a phase `theta` that turns at its natural frequency, `theta' = omega`.
Couplings to other oscillators add to the rate of the phase.
Args:
label:
The label of the oscillator or the group.
visual:
How the oscillator is drawn, or `None` for `default_visual_`.
theta:
The initial phase.
omega:
The natural frequency.
Attributes:
default_visual_:
The visual of the oscillator when it sets none.
Examples:
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.physics.common.oscillators import PhaseOscillatorElement
>>> system = System(elements=[PhaseOscillatorElement('clock', omega=2.0)])
>>> float(Simulation(system, dt=0.5).check().prepare().run(until=1.0).trajectory_.values('clock.theta')[-1])
2.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
theta: Annotated[float, StateSpec(role='phase')] = 0.0,
omega: Annotated[float, ParameterSpec()] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.theta = theta
self.omega = omega
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the turning of the phase at the natural frequency.
Args:
namespace:
The namespace of the symbols the laws can use.
Returns:
laws:
The rate of the phase.
"""
return [RateLaw(namespace.theta, namespace.omega)]
class KuramotoCouplingInteraction(CouplingInteraction):
"""Kuramoto coupling `theta_i' += K sum_j sin(theta_j - theta_i - alpha)` of the phases of a group.
The sum, or the mean, runs over the members `j` that the mode selects: the neighbours along the edges of a lattice
or network, or within a neighbourhood such as `RadiusNeighborhood(r)`, every other member of the group, every member
of a separate source group, or the next member. With `include_self` it also runs over `j = i`, whose term
`-K sin(alpha)` slows every phase. Over the whole group without weights or delay, the coupling reads the group sums
of `sin theta` and `cos theta` once, so it costs time linear in the size of the group.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the group of oscillators.
source:
The path of the group of oscillators the target couples to in the mode `'group'`, or `None` for the target.
K:
The coupling strength of one pair.
alpha:
The phase lag.
state:
The name of the state that holds the phases.
reduce:
`'sum'` or `'mean'`.
mode:
`'edges'`, `'group'`, or `'next'`.
weights:
The name of the state of the edge group in the mode `'edges'`, a matrix in the mode `'group'`, or `None`.
include_self:
Whether the sum also runs over the member itself.
delay:
The delay of the phases of the other members, 0 for none.
neighborhood:
The neighbourhood of the mode `'edges'`, or `None` for the edges of the lattice or network.
layer:
The layer of the network whose edges the interaction runs on, or `None` for the graph of the network.
Examples:
>>> import numpy as np
>>> from skcomplex.core.elements import GroupElement
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.physics.common.oscillators import (
... KuramotoCouplingInteraction,
... PhaseOscillatorElement,
... )
>>> phases = GroupElement('phases', element=PhaseOscillatorElement, size=4, values={'theta': [0, 1, 2, 3]})
>>> coupling = KuramotoCouplingInteraction('coupling', target='phases', K=2.0, mode='group', reduce='mean')
>>> system = System(elements=[phases], interactions=[coupling])
>>> theta = Simulation(system, dt=1.0).check().prepare().run(until=20.0).trajectory_.values('phases.theta')
>>> round(float(np.ptp(theta[-1])), 6)
0.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
source: Annotated[str | None, RoleSpec(optional=True)] = None,
K: Annotated[float, ParameterSpec()] = 1.0,
alpha: Annotated[float, ParameterSpec()] = 0.0,
state: str = 'theta',
reduce: str = 'sum',
mode: str = 'edges',
weights: str | npt.ArrayLike | scipy.sparse.sparray | None = None,
include_self: bool = False,
delay: float = 0.0,
neighborhood: BaseNeighborhood | None = None,
layer: str | None = None,
) -> None:
super().__init__(
label,
target=target,
source=source,
read=state,
write=state,
function=self._sine,
reduce=reduce,
mode=mode,
weights=weights,
include_self=include_self,
delay=delay,
neighborhood=neighborhood,
layer=layer,
)
self.K = K
self.alpha = alpha
self.state = state
def _sine(self, namespace: SimpleNamespace) -> Expression:
"""Return the term `K sin(theta_j - theta_i - alpha)` of one other member."""
member, other = getattr(namespace.member, self.state), getattr(namespace.other, self.state)
return namespace.K * sin(other - member - namespace.alpha)
def _get_group_reduction(self, namespace: SimpleNamespace) -> Target:
"""Return the reduction over the group, from the group sums of `sin theta` and `cos theta` when it can."""
if self.weights is not None or self.delay != 0:
return super()._get_group_reduction(namespace)
source, target = namespace.source, namespace.target
shifted = getattr(target, self.state) + namespace.alpha
phases = getattr(source, self.state)
sines, cosines = total(sin(phases), over=source), total(cos(phases), over=source)
reduction = namespace.K * (sines * cos(shifted) - cosines * sin(shifted))
if self._excludes_self():
reduction = reduction + namespace.K * sin(namespace.alpha)
if self.reduce == 'mean':
reduction = reduction / (count(over=source) - (1 if self._excludes_self() else 0))
return reduction
Equations
oscillators.theta[i]: diff(theta_oscillators(t), t) == (_reduced_1*cos(alpha_coupling + theta_oscillators(t)) - _reduced_2*sin(alpha_coupling + theta_oscillators(t)))*K_coupling/_reduced_0 + omega_oscillators