Skip to content

Lorenz-96

Lorenz-96, dx_i/dt = (x_{i+1} - x_{i-2}) x_{i-1} - x_i + F on a ring of sites around a latitude circle.

  • Name: lorenz_96
  • Field: physics/climate
  • Formalism: ordinary differential equations
  • References: Lorenz 1996; Datseris and Parlitz 2022

Factory

@register_model(
    field='physics/climate',
    tags=('chaos', 'ode', 'atmosphere', 'spatiotemporal chaos', 'data assimilation'),
    references=('Lorenz 1996', 'Datseris and Parlitz 2022'),
    defaults={'until': 20.0, 'dt': 0.05},
)
def lorenz_96(
    n_sites: int = 20,
    *,
    F: float = 8.0,
    perturbation: float = 0.01,
) -> System:
    """Return Lorenz-96, `dx_i/dt = (x_{i+1} - x_{i-2}) x_{i-1} - x_i + F` on a ring of sites around a latitude circle.

    The sites start at the steady state `F`, and the first one is displaced by `perturbation`. For `F = 8` the
    steady state is unstable and the ring is chaotic.

    Args:
        n_sites:
            The number of sites, at least 4, the dimension of the state `x`.

        F:
            The constant forcing.

        perturbation:
            The displacement of the first site from the steady state.

    Returns:
        system:
            The model, with the values `ring.x` of all sites.

    Examples:
        >>> from skcomplex.models.physics.climate import lorenz_96
        >>> lorenz_96(6).space.n_dim
        6
    """
    values = np.full(n_sites, float(F))
    values[0] += perturbation
    return System(elements=[Lorenz96Element('ring', x=values, F=F)], space=ContinuousSpace(n_dim=n_sites))

Components

class Lorenz96Element(Element):
    """Ring of Lorenz-96 sites around a latitude circle, advected by their neighbours, damped, and forced.

    Every site follows `dx_i/dt = (x_{i+1} - x_{i-2}) x_{i-1} - x_i + F`, with the indices taken around the ring.

    Args:
        label:
            The label of the ring.

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

        x:
            The initial values of all sites.

        F:
            The constant forcing.

    Examples:
        >>> from skcomplex.models.physics.climate import Lorenz96Element
        >>> Lorenz96Element('ring').F
        8.0
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        x: Annotated[npt.ArrayLike, StateSpec(vector=True)] = 8.0,
        F: Annotated[float, ParameterSpec()] = 8.0,
    ) -> None:
        super().__init__(label, visual=visual)
        self.x = x
        self.F = F

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the advection, damping, and forcing of every site of the ring.

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

        Returns:
            laws:
                The rate law of the values of all sites.
        """
        x = s.x
        n_sites = len(x)
        rates = [(x[(i + 1) % n_sites] - x[i - 2]) * x[i - 1] - x[i] + s.F for i in range(n_sites)]
        return [RateLaw(s.x, vector(SR, rates))]

Equations

ring.x[0]: diff(x_ring_0(t), t) == (x_ring_1(t) - x_ring_18(t))*x_ring_19(t) + F_ring - x_ring_0(t)
ring.x[1]: diff(x_ring_1(t), t) == -(x_ring_19(t) - x_ring_2(t))*x_ring_0(t) + F_ring - x_ring_1(t)
ring.x[2]: diff(x_ring_2(t), t) == -(x_ring_0(t) - x_ring_3(t))*x_ring_1(t) + F_ring - x_ring_2(t)
ring.x[3]: diff(x_ring_3(t), t) == -(x_ring_1(t) - x_ring_4(t))*x_ring_2(t) + F_ring - x_ring_3(t)
ring.x[4]: diff(x_ring_4(t), t) == -(x_ring_2(t) - x_ring_5(t))*x_ring_3(t) + F_ring - x_ring_4(t)
ring.x[5]: diff(x_ring_5(t), t) == -(x_ring_3(t) - x_ring_6(t))*x_ring_4(t) + F_ring - x_ring_5(t)
ring.x[6]: diff(x_ring_6(t), t) == -(x_ring_4(t) - x_ring_7(t))*x_ring_5(t) + F_ring - x_ring_6(t)
ring.x[7]: diff(x_ring_7(t), t) == -(x_ring_5(t) - x_ring_8(t))*x_ring_6(t) + F_ring - x_ring_7(t)
ring.x[8]: diff(x_ring_8(t), t) == -(x_ring_6(t) - x_ring_9(t))*x_ring_7(t) + F_ring - x_ring_8(t)
ring.x[9]: diff(x_ring_9(t), t) == (x_ring_10(t) - x_ring_7(t))*x_ring_8(t) + F_ring - x_ring_9(t)
ring.x[10]: diff(x_ring_10(t), t) == (x_ring_11(t) - x_ring_8(t))*x_ring_9(t) + F_ring - x_ring_10(t)
ring.x[11]: diff(x_ring_11(t), t) == (x_ring_12(t) - x_ring_9(t))*x_ring_10(t) + F_ring - x_ring_11(t)
ring.x[12]: diff(x_ring_12(t), t) == -(x_ring_10(t) - x_ring_13(t))*x_ring_11(t) + F_ring - x_ring_12(t)
ring.x[13]: diff(x_ring_13(t), t) == -(x_ring_11(t) - x_ring_14(t))*x_ring_12(t) + F_ring - x_ring_13(t)
ring.x[14]: diff(x_ring_14(t), t) == -(x_ring_12(t) - x_ring_15(t))*x_ring_13(t) + F_ring - x_ring_14(t)
ring.x[15]: diff(x_ring_15(t), t) == -(x_ring_13(t) - x_ring_16(t))*x_ring_14(t) + F_ring - x_ring_15(t)
ring.x[16]: diff(x_ring_16(t), t) == -(x_ring_14(t) - x_ring_17(t))*x_ring_15(t) + F_ring - x_ring_16(t)
ring.x[17]: diff(x_ring_17(t), t) == -(x_ring_15(t) - x_ring_18(t))*x_ring_16(t) + F_ring - x_ring_17(t)
ring.x[18]: diff(x_ring_18(t), t) == -(x_ring_16(t) - x_ring_19(t))*x_ring_17(t) + F_ring - x_ring_18(t)
ring.x[19]: diff(x_ring_19(t), t) == (x_ring_0(t) - x_ring_17(t))*x_ring_18(t) + F_ring - x_ring_19(t)