Van der Pol ice ages

Glacial cycles as a Van der Pol relaxation oscillator, optionally paced by orbital forcing.

• Name: van_der_pol_ice_ages
• Field: physics/climate
• Formalism: ordinary differential equations
• References: De Saedeleer, Crucifix, and Wieczorek 2013; Datseris and Parlitz 2022

Factory

@register_model(
    field='physics/climate',
    tags=('ode', 'limit cycle', 'relaxation oscillator', 'ice ages', 'paleoclimate'),
    references=('De Saedeleer, Crucifix, and Wieczorek 2013', 'Datseris and Parlitz 2022'),
    defaults={'until': 400.0, 'dt': 0.5},
)
def van_der_pol_ice_ages(
    x: float = 0.0,
    y: float = 1.0,
    *,
    alpha: float = 30.0,
    beta: float = 0.75,
    gamma: float = 0.0,
    tau: float = 35.09,
) -> System:
    """Return glacial cycles as a Van der Pol relaxation oscillator, optionally paced by orbital forcing.

    Time is in thousands of years. The slow variable `x` stands for ice volume and the fast variable `y` for the climate
    state, with `tau x' = -(y + beta + gamma F(t))` and `tau y' = -alpha (y^3/3 - y - x)`, where
    `F(t) = cos(2 pi t/23) + cos(2 pi t/41)` stands for precession and obliquity. For `|beta| < 1` the unforced
    oscillator has a relaxation limit cycle of slow ice build-up and fast terminations, and `gamma` sets how strongly
    the orbit paces it.

    Args:
        x:
            The initial ice volume.

        y:
            The initial climate state.

        alpha:
            The time-scale separation of the fast climate state.

        beta:
            The asymmetry of the cycle, between build-up and termination.

        gamma:
            The strength of the orbital forcing.

        tau:
            The time scale, in thousands of years.

    Returns:
        system:
            The model, with the states `climate.x` and `climate.y`.

    Examples:
        >>> from skcomplex.models.physics.climate import van_der_pol_ice_ages
        >>> van_der_pol_ice_ages(gamma=0.4).elements[0].gamma
        0.4
    """
    return System(
        elements=[VanDerPolIceAgesElement('climate', x=x, y=y, alpha=alpha, beta=beta, gamma=gamma, tau=tau)],
    )

Components

class VanDerPolIceAgesElement(Element):
    """Climate of glacial cycles as a Van der Pol relaxation oscillator, a slow ice volume and a fast climate state.

    With time in thousands of years, `tau x' = -(y + beta + gamma F(t))` and `tau y' = -alpha (y^3/3 - y - x)`, where
    `F(t) = cos(2 pi t/23) + cos(2 pi t/41)` stands for precession and obliquity.

    Args:
        label:
            The label of the climate.

        visual:
            How the element is drawn, or `None` for `default_visual_`.

        x:
            The initial ice volume.

        y:
            The initial climate state.

        alpha:
            The time-scale separation of the fast climate state.

        beta:
            The asymmetry of the cycle, between build-up and termination.

        gamma:
            The strength of the orbital forcing.

        tau:
            The time scale, in thousands of years.

    Examples:
        >>> from skcomplex.models.physics.climate import VanDerPolIceAgesElement
        >>> VanDerPolIceAgesElement('climate').alpha
        30.0
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        x: Annotated[float, StateSpec()] = 0.0,
        y: Annotated[float, StateSpec()] = 1.0,
        alpha: Annotated[float, ParameterSpec()] = 30.0,
        beta: Annotated[float, ParameterSpec()] = 0.75,
        gamma: Annotated[float, ParameterSpec()] = 0.0,
        tau: Annotated[float, ParameterSpec()] = 35.09,
    ) -> None:
        super().__init__(label, visual=visual)
        self.x = x
        self.y = y
        self.alpha = alpha
        self.beta = beta
        self.gamma = gamma
        self.tau = tau

    def _get_orbital_forcing(self, s: SimpleNamespace) -> Expression:
        """Return the orbital forcing, the sum of the cycles of precession and obliquity.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            forcing:
                The orbital forcing.
        """
        return cos(2 * pi * s.t / 23) + cos(2 * pi * s.t / 41)

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the slow ice volume and the fast climate of the relaxation oscillator under orbital forcing.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The rate laws of the ice volume and the climate state.
        """
        forcing = self._get_orbital_forcing(s)
        return [
            RateLaw(s.x, -(s.y + s.beta + s.gamma * forcing) / s.tau),
            RateLaw(s.y, -s.alpha * (s.y**3 / 3 - s.y - s.x) / s.tau),
        ]

Equations

climate.x: diff(x_climate(t), t) == -(gamma_climate*(cos(2/23*pi*t) + cos(2/41*pi*t)) + beta_climate + y_climate(t))/tau_climate
climate.y: diff(y_climate(t), t) == -1/3*(y_climate(t)^3 - 3*x_climate(t) - 3*y_climate(t))*alpha_climate/tau_climate