Lorenz 96 model

The Lorenz 96 model is a dynamical system formulated by Edward Lorenz in 1996.^[1] It is defined as follows. For $i=1,...,N$ :

{\frac {dx_{i}}{dt}}=(x_{i+1}-x_{i-2})x_{i-1}-x_{i}+F

where it is assumed that $x_{-1}=x_{N-1},x_{0}=x_{N}$ and $x_{N+1}=x_{1}$ and $N\geq 4$ . Here $x_{i}$ is the state of the system and $F$ is a forcing constant. $F=8$ is a common value known to cause chaotic behavior.

It is commonly used as a model problem in data assimilation.^[2]

Python simulation

from scipy.integrate import odeintimport matplotlib.pyplot as pltimport numpy as np# These are our constantsN = 5  # Number of variablesF = 8  # Forcingdef L96(x, t):    """Lorenz 96 model with constant forcing"""    # Setting up vector    d = np.zeros(N)    # Loops over indices (with operations and Python underflow indexing handling edge cases)    for i in range(N):        d[i] = (x[(i + 1) % N] - x[i - 2]) * x[i - 1] - x[i] + F    return dx0 = F * np.ones(N)  # Initial state (equilibrium)x0[0] += 0.01  # Add small perturbation to the first variablet = np.arange(0.0, 30.0, 0.01)x = odeint(L96, x0, t)# Plot the first three variablesfig = plt.figure()ax = fig.add_subplot(projection="3d")ax.plot(x[:, 0], x[:, 1], x[:, 2])ax.set_xlabel("$x_1$")ax.set_ylabel("$x_2$")ax.set_zlabel("$x_3$")plt.show()

Julia simulation

using DynamicalSystems, PyPlotPyPlot.using3D()# parameters and initial conditionsN = 5F = 8.0u₀ = F * ones(N)u₀[1] += 0.01 # small perturbation# The Lorenz-96 model is predefined in DynamicalSystems.jl:ds = Systems.lorenz96(N; F = F)# Equivalently, to define a fast version explicitly, do:struct Lorenz96{N} end # Structure for size typefunction (obj::Lorenz96{N})(dx, x, p, t) where {N}    F = p[1]    # 3 edge cases explicitly (performance)    @inbounds dx[1] = (x[2] - x[N - 1]) * x[N] - x[1] + F    @inbounds dx[2] = (x[3] - x[N]) * x[1] - x[2] + F    @inbounds dx[N] = (x[1] - x[N - 2]) * x[N - 1] - x[N] + F    # then the general case    for n in 3:(N - 1)      @inbounds dx[n] = (x[n + 1] - x[n - 2]) * x[n - 1] - x[n] + F    end    return nothingendlor96 = Lorenz96{N}() # create structds = ContinuousDynamicalSystem(lor96, u₀, [F])# And now evolve a trajectorydt = 0.01 # sampling timeTf = 30.0 # final timetr = trajectory(ds, Tf; dt = dt)# And plot in 3D:x, y, z = columns(tr)plot3D(x, y, z)

References

[1]

[2]