Lorenz 96 model

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

[math]\displaystyle{ \frac{dx_i}{dt} = (x_{i+1}-x_{i-2})x_{i-1} - x_i + F }[/math]

where it is assumed that [math]\displaystyle{ x_{-1}=x_{N-1},x_0=x_N }[/math] and [math]\displaystyle{ x_{N+1}=x_1 }[/math] and [math]\displaystyle{ N \ge 4 }[/math]. Here [math]\displaystyle{ x_i }[/math] is the state of the system and [math]\displaystyle{ F }[/math] is a forcing constant. [math]\displaystyle{ F=8 }[/math] is a common value known to cause chaotic behavior.

It is commonly used as a model problem in data assimilation.

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)