Lorenz 96 Model

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

:$\frac = (x_-x_)x_ - x_i + F$

where it is assumed that $x_=x_,x_0=x_N$ and $x_=x_1$ and $N \ge 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.

Python simulation

from mpl_toolkits.mplot3d import Axes3D
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 = (x i + 1) % N- x - 2 * x - 1- x + F
return d

x0 = F * np.ones(N) # Initial state (equilibrium)
x0 += 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.gca(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₀ += 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 end # Structure for size type
function (obj::Lorenz96)(dx, x, p, t) where
F = p # 3 edge cases explicitly (performance)
@inbounds dx = (x - x - 1 * x - x + F
@inbounds dx = (x - x * x - x + F
@inbounds dx = (x - x - 2 * x - 1- x + F
# then the general case
for n in 3:(N - 1)
@inbounds dx = (x + 1- x - 2 * x - 1- x + F
end
return nothing
end

lor96 = Lorenz96() # create struct
ds = ContinuousDynamicalSystem(lor96, u₀,

# 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

{{Reflist

Articles with example Python (programming language) code
Chaotic maps
Articles with example Julia code