# 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 {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}$ . 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 scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np

# These are our constants
N = 36  # Number of variables
F = 8  # Forcing

def Lorenz96(x, t):
# Compute state derivatives
d = np.zeros(N)
# First the 3 edge cases: i=1,2,N
d = (x - x[N-2]) * x[N-1] - x
d = (x - x[N-1]) * x - x
d[N-1] = (x - x[N-3]) * x[N-2] - x[N-1]
# Then the general case
for i in range(2, N-1):
d[i] = (x[i+1] - x[i-2]) * x[i-1] - x[i]
d = d + F

# Return the state derivatives
return d

x0 = F*np.ones(N)  # Initial state (equilibrium)
x0 += 0.01  # Add small perturbation to 20th variable
t = np.arange(0.0, 30.0, 0.01)

x = odeint(Lorenz96, x0, t)

# plot first three variables
from mpl_toolkits.mplot3d import Axes3D
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()