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

 

Plot of the first three variables of the simulation

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

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

1. Lorenz, Edward (1996). "Predictability – A problem partly solved" (PDF). Seminar on Predictability, Vol. I, ECMWF.
2. Ott, Edward; et al. (2002). "A Local Ensemble Kalman Filter for Atmospheric Data Assimilation". arXiv:physics/0203058.