The Lorenz system of coupled, ordinary, first-order differential equations have chaotic solutions for certain parameter values $\sigma$, $\rho$ and $\beta$ and initial conditions, $u(0)$, $v(0)$ and $w(0)$.

\begin{align*} \frac{\mathrm{d}u}{\mathrm{d}t} &= \sigma (v - u)\\ \frac{\mathrm{d}v}{\mathrm{d}t} &= \rho u - v - uw\\ \frac{\mathrm{d}w}{\mathrm{d}t} &= uv - \beta w \end{align*}

The following program plots the Lorenz attractor (the values of $x$, $y$ and $z$ as a parametric function of time) on a Matplotlib 3D projection.

*This code is also available on my github page.*

import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # Lorenz paramters and initial conditions sigma, beta, rho = 10, 2.667, 28 u0, v0, w0 = 0, 1, 1.05 # Maximum time point and total number of time points tmax, n = 100, 10000 def lorenz(X, t, sigma, beta, rho): """The Lorenz equations.""" u, v, w = X up = -sigma*(u - v) vp = rho*u - v - u*w wp = -beta*w + u*v return up, vp, wp # Integrate the Lorenz equations on the time grid t t = np.linspace(0, tmax, n) f = odeint(lorenz, (u0, v0, w0), t, args=(sigma, beta, rho)) x, y, z = f.T # Plot the Lorenz attractor using a Matplotlib 3D projection fig = plt.figure() ax = fig.gca(projection='3d') # Make the line multi-coloured by plotting it in segments of length s which # change in colour across the whole time series. s = 10 c = np.linspace(0,1,n) for i in range(0,n-s,s): ax.plot(x[i:i+s+1], y[i:i+s+1], z[i:i+s+1], color=(1,c[i],0), alpha=0.4) # Remove all the axis clutter, leaving just the curve. ax.set_axis_off() plt.show()

## Comments

Comments are pre-moderated. Please be patient and your comment will appear soon.

## Boris 2 years, 11 months ago

There is a discrepancy between the formula and the code for du/dt

Link | Reply## christian 2 years, 11 months ago

Thanks, Boris – there was a typo in the mark up for first equation of the Lorenz system. I've fixed it now.

Link | Reply## Vi 1 month, 2 weeks ago

Hello! Why 0, 1, 1.5 are the initial conditions?

Link | Reply## christian 1 month, 2 weeks ago

They don't have to be: it will attract from many different initial conditions (1, 1, 1) works, as does (1, 0, 0), etc.

Link | Reply## Asha. R 1 month, 2 weeks ago

Which method of solution in this program, for ex:RK4 like that which method?

Link | Reply## christian 1 month, 2 weeks ago

As written, the code uses scipy.optimize.odeint which wraps odepack's LSODA solver.

Link | Reply## New Comment