Almost every book on chaos or non-linear dynamics contain the Lorenz differential equations, prints the unmistakable image of the attractor, or mentions the "butterfly effect". It has been intensely studied by scientists from many disciplines since its publication by Edward Lorenz in 1963. In it, he described a simplified mathematical model for atmospheric convection [1] . Nowadays, these equations are simply known as Lorenz attractor. However, for many years scientist have argued if Lorenz attractor was truly chaos or an artifact of exponential and explosive amplifications of numerical truncation errors. The existence of Lorenz attractor was finally settled by Tucker in 2002 [2] .
Lorenz was a meteorologist and a mathematician in search of a model that was capable of generating long stretches of data statistically similar to atmospheric conditions. So, he set out to build a model with an aperiodic solution that was simple enough for the computing power available at the time. After a long search, Lorenz came up with the now-famous set of equations [3] :
where variables ($x$, $y$, $z$) and parameters ($\sigma$, $\rho$, $\beta$) are all real numbers.
Like many nonlinear equations, there isn't an analytical solution for Lorenz attractor. So, with practicality in mind, I will use numerical methods (more specifically discrete numerical methods with finite time-steps) to solve these equations and share some surprising features using the generated time-series. Besides, time-series are often the only observable data available to quantify the complicated dynamics in real life situations or experiments. So, let's jump right into it!
import utils
N=30
t, x_t = utils.solve_lorenz(10.0, 8.0/3, 28.0, N)
Above, we computed (numerically solved) Lorenz equation for $N=30$ trajectories with random initial conditions for $x$, $y$, $z$ using $\sigma = 10$, $\rho = 2.\overline{6}$, and $\beta = 28$. I will not get into what variables ($x$, $y$, $z$) and parameters ($\sigma$, $\rho$, $\beta$) represent. However, for wide range of values of parameters numerical solutions of the equations above generate extremely complicated trajectories. Let's plot them:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
# Pick a colormap: viridis, plasma, gray, binary, seismic, gnuplot
viridis = cm.get_cmap('viridis', N)
colors = viridis(np.linspace(0, 1, N))
plt.figure(figsize=(12, 12))
plt.style.use('bmh')
plt.xlabel('$x$')
plt.ylabel('$z$')
for i in range(N):
x, y, z = x_t[i,:,:].T
plt.plot(x, z, '-', c=colors[i])
This is the famous and unmistakable butterfly figure plotted in two dimensions. Notice that each trajectory exhibits aperiodic behavior, not cross itself if we consider the full three-dimensional picture.
Matplotlib is a powerful plotting library which produces publication quality figures. It can be very complex and time-consuming to create aesthetically pleasing plots. Luckily, we can take shortcuts and use ready-made styles. One particular style is called "bmh" which offer many improvements over the default settings of matplotlib and works well with Jupyter notebooks. This style comes from Bayesian Methods for Hackers [4] book. I find it particularly suiting for scientific graphing, and I will be using it frequently.
The aperiodic behavior is more apparent when we plot each variable separately.
plt.figure(figsize=(12, 18))
x, y, z = x_t[20,:,:].T
ax1 = plt.subplot(311)
ax1.plot(z)
ax1.set_ylabel('$z$')
ax2 = plt.subplot(312, sharex=ax1)
ax2.plot(y)
ax2.set_ylabel('$y$')
ax3 = plt.subplot(313, sharex=ax1)
ax3.plot(x)
ax3.set_ylabel('$x$')
ax3.set_xlabel('$n$')
plt.show()
The plots above show the variables ($x$, $y$, $z$) of a turbulent trajectory in three-dimensional space with the following properties:
For wide ranges of values of parameters ($\sigma$, $\rho$, $\beta$), the properties above hold true. The long time behavior of Lorenz equations result in a chaotic and unpredictable behavior. The trajectories are aperiodic and appears to loop around a " strange attractor " that resemble the wings of a butterfly. In other words, trajectories starting from two neighboring points, we can easily visualize and understand how the time evolution can be chaotic with extreme sensitivity to initial conditions.
Creating histograms of the average positions (across different trajectories) show that, on average, the trajectories swirl about the strange attractors.
t, x_t = utils.solve_lorenz(10.0, 8.0/3, 28.0, 100, 1000)
xyz_avg = x_t.mean(axis=1)
# xyz_avg.shape
plt.figure(figsize=(12, 4))
bx1 = plt.subplot(121)
bx1.hist(xyz_avg[:,0], 100)
bx1.set_xlabel('$x$')
bx1.set_ylabel('$n$')
bx1.set_title('Average $x$')
bx2 = plt.subplot(122)
bx2.hist(xyz_avg[:,1], 100)
bx2.set_xlabel('$z$')
bx1.set_ylabel('$n$')
bx2.set_title('Average $z$')
plt.show()
The histograms above show randomly selected 100 points ($x$, $y$, $z$) iterated for $t_{max}$ = 1000. If the trajectories were periodic, we would expect to see big spikes in the histograms representing limit cycles consisting of closed curves. Instead, we observe a broad distribution of points over a bounded range.
To better understand the intricate behavior of the Lorenz equation, we will use so-called Poincaré section in the next article .