The Lorenz equations¶
Here's a quick exploration of the Lorenz equations with Python. Of course, we'll use some libraries:
In [1]:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import root
from scipy.linalg import eig
Mathematically, the Lorenz equations are:
$$ \begin{align} x' &= \sigma(y-x) \\ y' &= x(\rho-z)-y \\ z' &= xy-\beta z. \end{align} $$The are many interesting values for the parameters; we'll take them to be
$$ \sigma = 10, \: \rho = 28, \text{ and } \: \beta=8/3. $$Here's a plot of a solution starting near the origin:
In [2]:
def lorenz(xyz, t=0, s=10, r=28, b=8/3):
x,y,z = xyz
xp = s*(y - x)
yp = r*x - y - x*z
zp = x*y - b*z
return xp, yp, zp
ts = np.linspace(0,100,10000)
sol = odeint(lorenz, [0,1,1.05], ts)
xs,ys,zs = sol.T
fig = plt.figure(figsize=(16,10))
ax = fig.gca(projection='3d')
ax.plot(xs, ys, zs, linewidth=0.3)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
plt.show()
We can plot this together with the equilibria:
In [3]:
fig = plt.figure(figsize=(16,10))
ax = fig.gca(projection='3d')
r0 = root(lorenz, [0,0,0]).x
r1 = root(lorenz, [5,5,20]).x
r2 = root(lorenz, [-5,-5,20]).x
ax.plot([r0[0],r1[0],r2[0]], [r0[1],r1[1],r2[1]],[r0[2],r1[2],r2[2]],'ko')
ax.plot(xs, ys, zs, linewidth=0.3)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
plt.show()
In order to understand what's going on near the equilibria, we should compute the eigenvalues of the Jacobian at those points:
In [4]:
def J(xyz, s=10, r=28, b=8/3):
x,y,z = xyz
return [
[-s, s, 0],
[r-z,-1,-x],
[y,x,-b]
]
e,v = eig(J(r1))
e
Out[4]:
There's your rotation! And here are two vectors determining the plane where the roation happens:
In [5]:
np.array([[z.real,z.imag] for z in v.transpose()[1]]).transpose()
Out[5]: