Kuramoto-Sakaguchi chimera
Two populations of identical phase oscillators with phase lag, where one synchronizes and one drifts.
- Name:
kuramoto_sakaguchi_chimera - Field:
physics/synchronization - Formalism: ordinary differential equations
- References: Abrams et al. 2008; Jensen 2023
Factory
@register_model(
field='physics/synchronization',
tags=('synchronization', 'oscillators', 'phase', 'chimera'),
references=('Abrams et al. 2008', 'Jensen 2023'),
defaults={'until': 100.0, 'dt': 1.0},
observables=lambda system: [*_order_parameter('first_order', 'first'), *_order_parameter('second_order', 'second')],
)
def kuramoto_sakaguchi_chimera(
n_oscillators: int = 32,
*,
A: float = 0.2,
beta: float = 0.1,
spread: float = 0.1,
seed: int = 0,
) -> System:
"""Return two populations of identical phase oscillators with phase lag, where one synchronizes and one drifts.
Oscillator `i` of population `p` follows `theta_i' = sum_q (K_pq / N) sum_j sin(theta_j - theta_i - alpha)`, with
the stronger coupling `mu = (1 + A) / 2` inside a population, the weaker `nu = (1 - A) / 2` between them, and the
lag `alpha = pi / 2 - beta`. Every population reads the means of `cos theta` and `sin theta` of both populations.
The first population starts nearly synchronized and the second spread over the circle. The populations settle in a
chimera, where one is fully synchronized and the other keeps drifting with a partial order parameter, which
breathes for `A = 0.2`. Which population ends up coherent depends on the initial phases.
Args:
n_oscillators:
The number of oscillators of every population.
A:
The disparity of the couplings inside and between the populations.
beta:
The distance of the phase lag from `pi / 2`.
spread:
The width of the initial phases of the first population.
seed:
The seed of the initial phases.
Returns:
system:
The model, with the phases `first.theta` and `second.theta` and the order parameters `first_order` and
`second_order`.
Examples:
>>> from skcomplex.models.physics.synchronization import kuramoto_sakaguchi_chimera
>>> len(kuramoto_sakaguchi_chimera(4).interactions)
4
"""
rng = np.random.default_rng(seed)
alpha, mu, nu = math.pi / 2 - beta, (1 + A) / 2, (1 - A) / 2
populations = [
GroupElement(
'first',
element=PhaseOscillatorElement,
size=n_oscillators,
values={'theta': rng.normal(0.0, spread, n_oscillators), 'omega': 0.0},
visual=VisualSpec('phase_ring', {'phase': 'theta'}),
),
GroupElement(
'second',
element=PhaseOscillatorElement,
size=n_oscillators,
values={'theta': _spread_phases(n_oscillators, seed + 1), 'omega': 0.0},
visual=VisualSpec('phase_ring', {'phase': 'theta'}),
),
]
interactions: list[BaseInteraction] = [
KuramotoCouplingInteraction(
f'{target}_from_{source}',
target=target,
source=source,
K=mu if target == source else nu,
alpha=alpha,
reduce='mean',
mode='group',
include_self=True,
)
for target in ('first', 'second')
for source in ('first', 'second')
]
return System(
elements=populations,
interactions=interactions,
)
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
first.theta[i]: diff(theta_first(t), t) == (_reduced_1*cos(alpha_first_from_first + theta_first(t)) - _reduced_2*sin(alpha_first_from_first + theta_first(t)))*K_first_from_first/_reduced_0 + (_reduced_4*cos(alpha_first_from_second + theta_first(t)) - _reduced_5*sin(alpha_first_from_second + theta_first(t)))*K_first_from_second/_reduced_3 + omega_first
second.theta[i]: diff(theta_second(t), t) == (_reduced_1*cos(alpha_second_from_first + theta_second(t)) - _reduced_2*sin(alpha_second_from_first + theta_second(t)))*K_second_from_first/_reduced_0 + (_reduced_4*cos(alpha_second_from_second + theta_second(t)) - _reduced_5*sin(alpha_second_from_second + theta_second(t)))*K_second_from_second/_reduced_3 + omega_second