Lorenz 96 model

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

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

where it is assumed that ${\displaystyle x_{-1}=x_{N-1},x_{0}=x_{N}}$ and ${\displaystyle x_{N+1}=x_{1}}$ and ${\displaystyle N\geq 4}$. Here ${\displaystyle x_{i}}$ is the state of the system and ${\displaystyle F}$ is a forcing constant. ${\displaystyle 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[edit]

from scipy.integrate import odeint import matplotlib.pyplot as plt import numpy as np  # These are our constants N = 5  # Number of variables F = 8  # Forcing   def 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 d   x0 = F * np.ones(N)  # Initial state (equilibrium) x0[0] += 0.01  # Add small perturbation to the first variable t = np.arange(0.0, 30.0, 0.01)  x = odeint(L96, x0, t)  # Plot the first three variables fig = 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[edit]

using DynamicalSystems, PyPlot PyPlot.using3D()  # parameters and initial conditions N = 5 F = 8.0 u₀ = 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 type function (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 nothing end  lor96 = Lorenz96{N}() # create struct ds = ContinuousDynamicalSystem(lor96, u₀, [F])  # And now evolve a trajectory dt = 0.01 # sampling time Tf = 30.0 # final time tr = trajectory(ds, Tf; dt = dt)  # And plot in 3D: x, y, z = columns(tr) plot3D(x, y, z) 

References[edit]