Winfree model

The Winfree model, theta_i' = omega_i + (K / N) sum_j P(theta_j) R(theta_i), with the sum over all j.

• Name: winfree
• Field: physics/synchronization
• Formalism: ordinary differential equations
• References: Winfree 1967; Ariaratnam and Strogatz 2001; Bick et al. 2020

Factory

@register_model(
    field='physics/synchronization',
    tags=('synchronization', 'oscillators', 'phase', 'oscillator death'),
    references=('Winfree 1967', 'Ariaratnam and Strogatz 2001', 'Bick et al. 2020'),
    defaults={'until': 50.0, 'dt': 0.5},
    observables=lambda system: _order_parameter('order', 'oscillators'),
)
def winfree(n_oscillators: int = 50, *, K: float = 0.5, spread: float = 0.1, seed: int = 0) -> System:
    """Return the Winfree model, `theta_i' = omega_i + (K / N) sum_j P(theta_j) R(theta_i)`, with the sum over all `j`.

    The influence `P(theta) = 1 + cos(theta)` and the sensitivity `R(theta) = -sin(theta)` make every oscillator respond
    to the mean influence of the population, which every oscillator reads as a group mean. The natural frequencies are
    spread evenly over `[1 - spread, 1 + spread]`. Weak coupling leaves the oscillators incoherent, moderate coupling
    locks them to one frequency, and strong coupling stops them, which is oscillator death.

    Args:
        n_oscillators:
            The number of oscillators.

        K:
            The coupling strength.

        spread:
            The half width of the interval of natural frequencies around 1.

        seed:
            The seed of the initial phases.

    Returns:
        system:
            The model, with the phases `oscillators.theta` and the order parameter `order`.

    Examples:
        >>> from skcomplex.models.physics.synchronization import winfree
        >>> len(winfree(5).interactions)
        1
    """
    oscillators = GroupElement(
        'oscillators',
        element=PhaseOscillatorElement,
        size=n_oscillators,
        values={
            'theta': np.random.default_rng(seed).uniform(0, 2 * np.pi, n_oscillators),
            'omega': np.linspace(1 - spread, 1 + spread, n_oscillators),
        },
    )
    return System(
        elements=[oscillators],
        interactions=[WinfreeMeanFieldInteraction('coupling', target='oscillators', K=K)],
    )

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 WinfreeMeanFieldInteraction(CouplingInteraction):
    """Winfree coupling `theta_i' += K mean_j P(theta_j) R(theta_i)` over every member `j` of the group.

    The influence `P(theta) = 1 + cos(theta)` that an oscillator sends and the sensitivity `R(theta) = -sin(theta)` with
    which it responds make every oscillator respond to the mean influence 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.

        K:
            The coupling strength.

    Examples:
        >>> from types import SimpleNamespace
        >>> from sage.all__sagemath_symbolics import var
        >>> from skcomplex.models.physics.synchronization import WinfreeMeanFieldInteraction
        >>> coupling = WinfreeMeanFieldInteraction('coupling', target='oscillators')
        >>> member, other = SimpleNamespace(theta=var('a')), SimpleNamespace(theta=var('b'))
        >>> coupling.function(SimpleNamespace(K=var('K'), member=member, other=other))
        -K*(cos(b) + 1)*sin(a)
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        K: Annotated[float, ParameterSpec()] = 1.0,
    ) -> None:
        super().__init__(
            label,
            target=target,
            read='theta',
            write='theta',
            function=self._response,
            reduce='mean',
            mode='group',
            include_self=True,
        )
        self.K = K

    def _response(self, namespace: SimpleNamespace) -> Expression:
        """Return the response `K P(theta_j) R(theta_i)` of a member to one other member."""
        return namespace.K * (1 + cos(namespace.other.theta)) * -sin(namespace.member.theta)

Equations

oscillators.theta[i]: diff(theta_oscillators(t), t) == omega_oscillators - (K_coupling*_reduced_0*sin(theta_oscillators(t)) + K_coupling*_reduced_1*sin(theta_oscillators(t)))/_reduced_2