Globally coupled Stuart-Landau oscillators with frequency spread
Globally coupled Stuart-Landau oscillators, whose oscillations all die for strong coupling and spread.
- Name:
oscillator_death - Field:
physics/synchronization - Formalism: ordinary differential equations
- References: Ermentrout 1990; Mirollo and Strogatz 1990; Bertin 2021
Factory
@register_model(
field='physics/synchronization',
tags=('synchronization', 'oscillators', 'amplitude death', 'mean field', 'ode'),
references=('Ermentrout 1990', 'Mirollo and Strogatz 1990', 'Bertin 2021'),
defaults={'until': 20.0, 'dt': 0.2},
observables=lambda system: [
ObservableSpec('amplitude', _amplitude, group='oscillators'),
ObservableSpec('X', group='oscillators', variable='x'),
ObservableSpec('Y', group='oscillators', variable='y'),
ObservableSpec('order', _order),
],
)
def oscillator_death(
n_oscillators: int = 100,
*,
K: float = 3.0,
width: float = 1.5,
distribution: str = 'lorentzian',
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return globally coupled Stuart-Landau oscillators, whose oscillations all die for strong coupling and spread.
Oscillator `j` follows `z_j' = (1 + iω_j) z_j - |z_j|² z_j + (K / N) Σ_k (z_k - z_j)` for `z_j = x_j + i y_j`.
Alone it settles on the unit circle. The frequencies are the midpoint quantiles of a density of width `ω0`, the
Lorentzian `ω0 / (π (ω² + ω0²))` or the quartic `g0 / (ω⁴ + ω0⁴)`. For many oscillators, the common fixed point
`z = 0` is stable for `K > 1` and `ω0` above a critical width, 1 for the Lorentzian and
`(1 + √(2K - 1)) / √2` for the quartic density, and every amplitude then decays to zero. For the Lorentzian, the
slowest perturbation decays at the rate `ω0 - 1`.
Args:
n_oscillators:
The number of oscillators.
K:
The coupling strength.
width:
The width `ω0` of the density of frequencies.
distribution:
The density of frequencies, `'lorentzian'` or `'quartic'`.
random_state:
The seed or generator of the initial phases.
Returns:
system:
The model, with the states `oscillators.x` and `oscillators.y` and the observables `amplitude`, the mean
amplitude, `X` and `Y`, the mean field, and `order`, its modulus.
Raises:
ValidationError:
The distribution is neither `'lorentzian'` nor `'quartic'`.
Examples:
>>> from skcomplex.models.physics.synchronization import oscillator_death
>>> oscillator_death(4, distribution='quartic').elements[0].check().prepare().size_
4
"""
omegas = _frequency_quantiles(n_oscillators, width, distribution)
x, y = _spread_on_circle(n_oscillators, random_state)
oscillators = GroupElement(
'oscillators',
element=StuartLandauElement,
size=n_oscillators,
values={'x': x, 'y': y, 'mu': 1.0, 'omega': omegas, 'tau': 1.0},
)
return System(
elements=[oscillators],
interactions=[
GlobalCouplingInteraction(f'coupling_{state}', target='oscillators', state=state, K=K) for state in 'xy'
],
)
Components
class StuartLandauElement(Element):
"""Stuart-Landau oscillator `z' = tau ((mu + i omega) z - kappa (1 + i beta) |z|^2 z)` of `z = x + i y`.
It is the normal form of a supercritical Hopf bifurcation. For `mu > 0` the oscillator settles on the limit cycle
of radius `sqrt(mu / kappa)`, which it turns at the frequency `tau (omega - beta mu)`, and for `mu < 0` it decays to
the fixed point 0. The twist `beta` makes the frequency depend on the amplitude, and the time scale `tau` speeds
up or slows down the whole flow.
Args:
label:
The label of the oscillator or the group.
visual:
How the oscillator is drawn, or `None` for `default_visual_`.
x:
The initial real part of `z`.
y:
The initial imaginary part of `z`.
mu:
The growth rate of the fixed point 0, the bifurcation parameter.
omega:
The frequency of the linear rotation.
beta:
The twist, the shear of the frequency with the amplitude.
kappa:
The strength of the cubic saturation.
tau:
The time scale.
Attributes:
default_visual_:
The visual of the oscillator when it sets none.
Examples:
>>> import math
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.physics.common.oscillators import StuartLandauElement
>>> system = System(elements=[StuartLandauElement('oscillator', x=0.1, mu=4.0)])
>>> trajectory = Simulation(system, dt=0.5).check().prepare().run(until=10.0).trajectory_
>>> round(math.hypot(trajectory.values('oscillator.x')[-1], trajectory.values('oscillator.y')[-1]), 4)
2.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
x: Annotated[float, StateSpec()] = 1.0,
y: Annotated[float, StateSpec()] = 0.0,
mu: Annotated[float, ParameterSpec()] = 1.0,
omega: Annotated[float, ParameterSpec()] = 1.0,
beta: Annotated[float, ParameterSpec()] = 0.0,
kappa: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.x = x
self.y = y
self.mu = mu
self.omega = omega
self.beta = beta
self.kappa = kappa
self.tau = tau
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the rates of the real and imaginary parts of `z`.
Args:
namespace:
The namespace of the symbols the laws can use.
Returns:
laws:
The rates of `x` and `y`.
"""
s = namespace
saturation = s.kappa * (s.x**2 + s.y**2)
return [
RateLaw(s.x, s.tau * (s.mu * s.x - s.omega * s.y - saturation * (s.x - s.beta * s.y))),
RateLaw(s.y, s.tau * (s.mu * s.y + s.omega * s.x - saturation * (s.y + s.beta * s.x))),
]
class GlobalCouplingInteraction(CouplingInteraction):
"""Diffusive coupling of one state of every member to the mean of the group, `x_i' += K (mean_j x_j - x_i)`.
The mean runs over every member, the member itself included.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the group.
state:
The name of the coupled state.
K:
The coupling strength.
Examples:
>>> from types import SimpleNamespace
>>> from sage.all__sagemath_symbolics import var
>>> from skcomplex.models.physics.synchronization import GlobalCouplingInteraction
>>> coupling = GlobalCouplingInteraction('coupling', target='oscillators', state='y').check().prepare()
>>> member, other = SimpleNamespace(y=var('y')), SimpleNamespace(y=var('z'))
>>> coupling.function(SimpleNamespace(K=var('K'), member=member, other=other))
-K*(y - z)
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'x',
K: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
super().__init__(
label,
target=target,
read=state,
write=state,
function=self._pull,
reduce='mean',
mode='group',
include_self=True,
)
self.state = state
self.K = K
def _pull(self, namespace: SimpleNamespace) -> Expression:
"""Return the pull `K (x_j - x_i)` of one other member."""
return namespace.K * (getattr(namespace.other, self.state) - getattr(namespace.member, self.state))
Equations
oscillators.x[i]: diff(x_oscillators(t), t) == ((x_oscillators(t)^2 + y_oscillators(t)^2)*(beta_oscillators*y_oscillators(t) - x_oscillators(t))*kappa_oscillators + mu_oscillators*x_oscillators(t) - omega_oscillators*y_oscillators(t))*tau_oscillators - (K_coupling_x*_reduced_0*x_oscillators(t) - K_coupling_x*_reduced_1)/_reduced_2
oscillators.y[i]: diff(y_oscillators(t), t) == -((beta_oscillators*x_oscillators(t) + y_oscillators(t))*(x_oscillators(t)^2 + y_oscillators(t)^2)*kappa_oscillators - omega_oscillators*x_oscillators(t) - mu_oscillators*y_oscillators(t))*tau_oscillators - (K_coupling_y*_reduced_0*y_oscillators(t) - K_coupling_y*_reduced_3)/_reduced_2