Globally coupled systems with heterogeneous time scales
Globally coupled oscillators with heterogeneous time scales, whose fixed point the coupling restabilizes.
- Name:
heterogeneous_global_coupling - Field:
physics/synchronization - Formalism: ordinary differential equations
- References: De Monte, d'Ovidio, and Mosekilde 2003; Bertin 2021
Factory
@register_model(
field='physics/synchronization',
tags=('synchronization', 'oscillators', 'amplitude death', 'mean field', 'heterogeneity', 'ode'),
references=("De Monte, d'Ovidio, and Mosekilde 2003", 'Bertin 2021'),
defaults={'until': 40.0, 'dt': 0.2},
observables=lambda system: [
ObservableSpec('amplitude', _amplitude, group='systems'),
ObservableSpec('X', group='systems', variable='x'),
ObservableSpec('Y', group='systems', variable='y'),
ObservableSpec('V', _time_scale_deviation_x, group='systems'),
ObservableSpec('W', _time_scale_deviation_y, group='systems'),
],
)
def heterogeneous_global_coupling(
n_systems: int = 100,
*,
k: float = 2.0,
beta: float = 0.1,
gamma: float = 2.0,
heterogeneity: float = 0.8,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return globally coupled oscillators with heterogeneous time scales, whose fixed point the coupling restabilizes.
System `j` follows `x_j' = τ_j g(x_j, y_j) + k (X - x_j)` and `y_j' = τ_j h(x_j, y_j) + k (Y - y_j)`, where `X`
and `Y` are the population means and `g + ih = (beta + i gamma) z - |z|² z` is the Hopf normal form of
`z = x + iy`. Its fixed point 0 is an unstable focus with eigenvalues `beta ± i gamma`. The time scales `τ_j` are
spread evenly over `[1 - heterogeneity, 1 + heterogeneity]`, with mean 1 and variance `sigma²`. The order
parameters `X`, `Y`, `V = <δτ δx>`, and `W = <δτ δy>` close on four equations. The global fixed point is stable
when their eigenvalues `nu - k / 2 ± √(k² / 4 + sigma² nu²)`, for `nu = beta ± i gamma`, have negative real parts,
which holds roughly for `2 beta (1 + sigma) < k < (gamma² - beta²) sigma² / beta`.
Args:
n_systems:
The number of coupled systems.
k:
The coupling constant.
beta:
The growth rate `β` of the unstable focus.
gamma:
The rotation rate `gamma` of the unstable focus.
heterogeneity:
The half width of the interval of time scales, below 1.
random_state:
The seed or generator of the initial phases.
Returns:
system:
The model, with the states `systems.x` and `systems.y` and the observables `amplitude`, the mean amplitude,
and the order parameters `X`, `Y`, `V`, and `W`.
Examples:
>>> from skcomplex.models.physics.synchronization import heterogeneous_global_coupling
>>> heterogeneous_global_coupling(4).elements[0].values['tau']
array([0.4, 0.8, 1.2, 1.6])
"""
taus = 1 + heterogeneity * (2 * (np.arange(n_systems) + 1 / 2) / n_systems - 1)
x, y = _spread_on_circle(n_systems, random_state)
radius = math.sqrt(abs(beta))
systems = GroupElement(
'systems',
element=StuartLandauElement,
size=n_systems,
values={'x': radius * x, 'y': radius * y, 'mu': beta, 'omega': gamma, 'tau': taus},
)
return System(
elements=[systems],
interactions=[
GlobalCouplingInteraction(f'coupling_{state}', target='systems', 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
systems.x[i]: diff(x_systems(t), t) == ((x_systems(t)^2 + y_systems(t)^2)*(beta_systems*y_systems(t) - x_systems(t))*kappa_systems + mu_systems*x_systems(t) - omega_systems*y_systems(t))*tau_systems - (K_coupling_x*_reduced_0*x_systems(t) - K_coupling_x*_reduced_1)/_reduced_2
systems.y[i]: diff(y_systems(t), t) == -((beta_systems*x_systems(t) + y_systems(t))*(x_systems(t)^2 + y_systems(t)^2)*kappa_systems - omega_systems*x_systems(t) - mu_systems*y_systems(t))*tau_systems - (K_coupling_y*_reduced_0*y_systems(t) - K_coupling_y*_reduced_3)/_reduced_2