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 solve_ivp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Create an image of the Lorenz attractor.
# The maths behind this code is described in the scipython blog article
# at https://scipython.com/blog/the-lorenz-attractor/
# Christian Hill, January 2016.
# Updated, January 2021 to use scipy.integrate.solve_ivp.
WIDTH, HEIGHT, DPI = 1000, 750, 100
# 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(t, X, 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.
soln = solve_ivp(lorenz, (0, tmax), (u0, v0, w0), args=(sigma, beta, rho),
dense_output=True)
# Interpolate solution onto the time grid, t.
t = np.linspace(0, tmax, n)
x, y, z = soln.sol(t)
# Plot the Lorenz attractor using a Matplotlib 3D projection.
fig = plt.figure(facecolor='k', figsize=(WIDTH/DPI, HEIGHT/DPI))
ax = fig.gca(projection='3d')
ax.set_facecolor('k')
fig.subplots_adjust(left=0, right=1, bottom=0, top=1)
# Make the line multi-coloured by plotting it in segments of length s which
# change in colour across the whole time series.
s = 10
cmap = plt.cm.winter
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=cmap(i/n), alpha=0.4)
# Remove all the axis clutter, leaving just the curve.
ax.set_axis_off()
plt.savefig('lorenz.png', dpi=DPI)
plt.show()
```

## Comments

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

## Boris 5 years, 11 months ago

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

Link | Reply## christian 5 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 3 years ago

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

Link | Reply## christian 3 years 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## Dr Dre 1 year, 6 months ago

Is there are reason why these initial conditions work?

Link | Reply## christian 1 year, 6 months ago

These are just the initial conditions that Lorenz himself chose. You get slightly different shapes with nearby values, and different behaviour (e.g. non-chaotic behaviour) with values far away. You can read about this on Wikipedia: https://en.wikipedia.org/wiki/Lorenz_system

Link | Reply## Asha. R 3 years ago

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

Link | Reply## christian 3 years ago

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

Link | Reply## Starslayerx 2 years, 4 months ago

why do n-s in range(0,n-s,s):

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

This is just for colouring the line: divide it into n/s segments and colour points 0->s one colour, s->2s the next colour, and so on until the points n-s -> n. I think the "+1" shouldn't be there in the plot command and is the reason for the dots: these are points that are plotted twice.

Link | Reply## New Comment