Delayed Kuramoto

Kuramoto phase oscillators coupled all to all through signals that arrive after a delay.

• Name: delayed_kuramoto
• Field: physics/synchronization
• Formalism: ordinary differential equations
• References: Yeung and Strogatz 1999; Kuramoto 1984

Factory

@register_model(
    field='physics/synchronization',
    tags=('synchronization', 'delay', 'oscillators'),
    references=('Yeung and Strogatz 1999', 'Kuramoto 1984'),
    defaults={'until': 50.0, 'dt': 0.1},
)
def delayed_kuramoto(
    n_oscillators: int = 5,
    *,
    K: float = 2.0,
    tau: float = 1.0,
    spread: float = 0.2,
    seed: int = 0,
) -> System:
    """Return Kuramoto phase oscillators coupled all to all through signals that arrive after a delay.

    Oscillator `i` follows `theta_i' = omega_i + (K / N) sum_j sin(theta_j(t - tau) - theta_i)`, with natural
    frequencies spread evenly around 1 and initial phases spread around the circle.

    Args:
        n_oscillators:
            The number of oscillators.

        K:
            The coupling strength.

        tau:
            The delay of the coupling.

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

        seed:
            The seed of the initial phases.

    Returns:
        system:
            The model, with the phases `oscillator_<i>.theta`.

    Examples:
        >>> from skcomplex.models.physics.synchronization import delayed_kuramoto
        >>> len(delayed_kuramoto(n_oscillators=3).interactions)
        6
    """
    phases = [2 * math.pi * ((seed * 0.618 + index / n_oscillators) % 1.0) for index in range(n_oscillators)]
    omegas = [1.0 + spread * (index / max(n_oscillators - 1, 1) - 0.5) for index in range(n_oscillators)]
    elements = [
        PhaseOscillatorElement(f'oscillator_{index}', theta=phases[index], omega=omegas[index])
        for index in range(n_oscillators)
    ]
    interactions = [
        DelayedCouplingInteraction(
            f'coupling_{source}_{target}',
            source=f'oscillator_{source}',
            target=f'oscillator_{target}',
            K=K / n_oscillators,
            tau=tau,
        )
        for source in range(n_oscillators)
        for target in range(n_oscillators)
        if source != target
    ]
    return System(elements=elements, 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 DelayedCouplingInteraction(RateInteraction):
    """Delayed sine coupling of one phase oscillator to another, `theta_t' += K sin(theta_s(t - tau) - theta_t)`.

    The target reads the phase of the source as it was a delay `tau` earlier.

    Args:
        label:
            The label of the interaction, unique among the parts of its system.

        source:
            The path of the oscillator whose delayed phase the coupling reads.

        target:
            The path of the oscillator whose phase the coupling turns.

        K:
            The coupling strength.

        tau:
            The delay.

    Examples:
        >>> from types import SimpleNamespace
        >>> from sage.all__sagemath_symbolics import var
        >>> from skcomplex.models.physics.synchronization import DelayedCouplingInteraction
        >>> coupling = DelayedCouplingInteraction('coupling', source='a', target='b').check().prepare()
        >>> namespace = SimpleNamespace(source=SimpleNamespace(theta=var('a')), target=SimpleNamespace(theta=var('b')))
        >>> [law.value for law in coupling.run(SimpleNamespace(**vars(namespace), K=var('K'), tau=1)).laws_]
        [K*sin(-b + delayed(a, 1))]
    """

    def __init__(
        self,
        label: str,
        *,
        source: Annotated[str, RoleSpec()],
        target: Annotated[str, RoleSpec()],
        K: Annotated[float, ParameterSpec()] = 1.0,
        tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    ) -> None:
        super().__init__(label, target=target, state='theta')
        self.source = source
        self.K = K
        self.tau = tau

    def _get_rate(self, namespace: SimpleNamespace) -> Expression:
        """Return `K sin(theta_s(t - tau) - theta_t)`.

        Args:
            namespace:
                The namespace of the symbols the rate can use.

        Returns:
            rate:
                The rate of the phase of the target.
        """
        return namespace.K * sin(delayed(namespace.source.theta, namespace.tau) - namespace.target.theta)

Equations

oscillator_0.theta: diff(theta_oscillator_0(t), t) == K_coupling_1_0*sin(-theta_oscillator_0(t) + theta_oscillator_1(t - tau_coupling_1_0)) + K_coupling_2_0*sin(-theta_oscillator_0(t) + theta_oscillator_2(t - tau_coupling_2_0)) + K_coupling_3_0*sin(-theta_oscillator_0(t) + theta_oscillator_3(t - tau_coupling_3_0)) + K_coupling_4_0*sin(-theta_oscillator_0(t) + theta_oscillator_4(t - tau_coupling_4_0)) + omega_oscillator_0
oscillator_1.theta: diff(theta_oscillator_1(t), t) == K_coupling_0_1*sin(theta_oscillator_0(t - tau_coupling_0_1) - theta_oscillator_1(t)) + K_coupling_2_1*sin(-theta_oscillator_1(t) + theta_oscillator_2(t - tau_coupling_2_1)) + K_coupling_3_1*sin(-theta_oscillator_1(t) + theta_oscillator_3(t - tau_coupling_3_1)) + K_coupling_4_1*sin(-theta_oscillator_1(t) + theta_oscillator_4(t - tau_coupling_4_1)) + omega_oscillator_1
oscillator_2.theta: diff(theta_oscillator_2(t), t) == K_coupling_0_2*sin(theta_oscillator_0(t - tau_coupling_0_2) - theta_oscillator_2(t)) + K_coupling_1_2*sin(theta_oscillator_1(t - tau_coupling_1_2) - theta_oscillator_2(t)) + K_coupling_3_2*sin(-theta_oscillator_2(t) + theta_oscillator_3(t - tau_coupling_3_2)) + K_coupling_4_2*sin(-theta_oscillator_2(t) + theta_oscillator_4(t - tau_coupling_4_2)) + omega_oscillator_2
oscillator_3.theta: diff(theta_oscillator_3(t), t) == K_coupling_0_3*sin(theta_oscillator_0(t - tau_coupling_0_3) - theta_oscillator_3(t)) + K_coupling_1_3*sin(theta_oscillator_1(t - tau_coupling_1_3) - theta_oscillator_3(t)) + K_coupling_2_3*sin(theta_oscillator_2(t - tau_coupling_2_3) - theta_oscillator_3(t)) + K_coupling_4_3*sin(-theta_oscillator_3(t) + theta_oscillator_4(t - tau_coupling_4_3)) + omega_oscillator_3
oscillator_4.theta: diff(theta_oscillator_4(t), t) == K_coupling_0_4*sin(theta_oscillator_0(t - tau_coupling_0_4) - theta_oscillator_4(t)) + K_coupling_1_4*sin(theta_oscillator_1(t - tau_coupling_1_4) - theta_oscillator_4(t)) + K_coupling_2_4*sin(theta_oscillator_2(t - tau_coupling_2_4) - theta_oscillator_4(t)) + K_coupling_3_4*sin(theta_oscillator_3(t - tau_coupling_3_4) - theta_oscillator_4(t)) + omega_oscillator_4