The Butterfly Effect: Meet the Lorenz Attractor
The Butterfly effect is more often than not misunderstood as the adage that the flap of a butterfly’s wings can cause a hurricane. While this is accurate in the macroscopic sense, it is certainly not the full picture. Here’s how the story goes.
Enter the Lorenz equations. The Lorenz equations were developed by Edward Lorenz in an attempt to model atmospheric convection — a phenomenon that deals with the movement of fluid parcels as a consequence of temperature fluctuations. The system consists of three coupled, first-order, nonlinear differential equations that are the hallmark of chaos. Lorenz described chaos as “when the present determines the future, but the approximate present does not approximately determine the future.” He found that nearly indistinguishable initial conditions could produce completely divergent outcomes, rendering weather prediction impossible beyond a time horizon of about a fortnight.
The Lorenz attractor is a set of chaotic solutions of the Lorenz system and is possibly the most famous depiction of chaotic behavior (save for the double pendulum, of course). The system itself describes the movement of a point in a three-dimensional space over time; it is formally described by the spatial coordinates x, y, and z, along with the constant system parameters of sigma, rho, and beta.
For his famous depiction of chaos, Lorenz used the values sigma = 10, beta = 8/3 and rho = 28. In the neighborhood of these parameter values, the system exhibits deterministic chaos. These trajectories never overlap and the system never lands on the same point twice, due to its fractal geometry. Now, plain old attractors are a set of states toward which a system tends to evolve, for a wide variety of starting conditions of the system. This translates to achieving approximately the same result as long as the starting conditions are approximately the same as well — some initial conditions x, y, and z should yield the same final result as some initial conditions x + e, y + e, and z + e where e is some small quantity. With strange attractors, however, this is not the case—the presence of the slightest bit of noise results in amplification in diverging trajectories as time progresses. The system is inherently unstable — approximate initial conditions do not approximately yield the same result.
Well, now, naturally one is forced to wonder whether trajectories go out to infinity. If the system is indeed inherently unstable, would trajectories eventually explode? The answer’s no. We don’t need the math for this. We can simply compute the vector field on a sphere and it suffices to see that everything points inwards. How to do this? Exercise left to the reader. You’re welcome.
Let’s get to the good stuff now. The general Lorenz attractor is covered in the matplotlib documentation. We aren’t interested in solely this. No. We’re better than that. We’ll begin by using Pillow to convert a series of frames into a .gif suitable for demonstrating the evolution of the Lorenz attractor. In the next article, we will also demonstrate this in Manim because the animation would be beautiful. I pity my laptop for the abuse I put it through in rendering this Manimation. The current render time is in the neighborhood of approximately two hours. Anyway, here’s to Python first.
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt, IPython.display as IPdisplay
import glob, os
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
from PIL import Image
save_folder = 'images/lorenz' # folder destination
if not os.path.exists(save_folder):
os.makedirs(save_folder)
# initial values (the starting x, y, z positions in space).
initial_state = [0.5, 0, 0]
# define the system parameters sigma, rho, and beta
sigma = 10
rho = 28
beta = 8/3
# define the time points to solve for
start_time = 1
end_time = 80
interval = 200 #FIX ME! a large interval will yield a smoother .gif
time_points = np.linspace(start_time, end_time, end_time * interval)
# define the lorenz system
def lorenz_system(current_state, t):
current_state = x, y, z
xdot = sigma * (y - x)
ydot = x * (rho - z) - y
zdot = x * y - beta * z
return [xdot, ydot, zdot] #returns the time derivatives
# plot the system in 3 dimensions
def plot_lorenz(xyz, n, save_folder, title_font):
fig = plt.figure(figsize=(12, 9))
x = xyz[:, 0]
y = xyz[:, 1]
z = xyz[:, 2]
ax.plot(x, y, z, color='purple', alpha=0.7, linewidth=0.7)
ax.set_xlim((-30, 30))
ax.set_ylim((-30, 30))
ax.set_zlim((0, 50))
plt.savefig(f'{save_folder}/{n:03d}.png', dpi=60,
bbox_inches='tight', pad_inches=0.1)
plt.close()
# return a list in iteratively larger discrete units
def get_units(list, size):
size = max(1, size)
chunks = [full_list[0:i] for i in range(1, len(full_list) + 1, size)]
return units
# get incrementally larger time units to reveal the attractor one frame at a time
units = get_units(time_points, size=20)
points = [odeint(lorenz_system, initial_state, unit) for unit in units]
#plotting points
for n, point in enumerate(points):
plot_lorenz(point, n, save_folder)
# create a tuple of display durations, one for each frame
firstlast = 100 #show the first and last frames for 100 ms
standard_duration = 5 #show all other frames for 5 ms
#subtract two for the first/last frames
durations = tuple([first_last] + [standard_duration] * (len(points) - 2) + [firstlast])
images = [Image.open(image) for image in glob.glob('{}/*.png'.format(save_folder))]
gif_filepath = 'images/animated-lorenz-attractor.gif'
# save as an animated gif
gif = images[0]
gif.info['duration'] = durations #ms per frame
gif.info['loop'] = 0 #how many times to loop (0=infinite)
gif.save(fp=gif_filepath, format='gif', save_all=True, append_images=images[1:])
IPdisplay.Image(url=gif_filepath)Below is the final result.
Let’s now compute the cross-sectional profiles.
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# the same run from the previous example
initial_state = [0, 1, 0]
sigma = 10
rho = 28
beta = 8 / 3
# define the time points to solve for, evenly spaced between the start and end times
start_time = 0
end_time = 100
time_points = np.linspace(start_time, end_time, end_time * 100)
# define the lorenz system
def lorenz_system(current_state, t):
x, y, z = current_state
xdot = sigma * (y - x)
ydot = x * (rho - z) - y
zdot = x * y - beta * z
return [xdot, ydot, zdot] # returns the time derivatives
# use odeint() to solve a system of ordinary differential equations
# the arguments are:
# 1, a function - computes the derivatives
# 2, a vector of initial system conditions (aka x, y, z positions in space)
# 3, a sequence of time points to solve for
# returns an array of x, y, and z value arrays for each time point, with the initial values in the first row
xyz = odeint(lorenz_system, initial_state, time_points)
# extract the individual arrays of x, y, and z values from the array of arrays
x = xyz[:, 0]
y = xyz[:, 1]
z = xyz[:, 2]
# plot the lorenz attractor in three-dimensional phase space
fig = plt.figure(figsize=(12, 9))
ax = fig.add_subplot(111, projection='3d')
ax.set_xlim((-30, 30))
ax.set_ylim((-30, 30))
ax.set_zlim((0, 50))
ax.plot(x, y, z, color='purple', alpha=0.7, linewidth=0.3)
ax.set_title('Lorenz attractor phase diagram')
plt.show()
For the two-dimensional cuts, we use the following.
# now plot two-dimensional cuts of the three-dimensional phase space
fig, ax = plt.subplots(1, 3, sharex=False, sharey=False, figsize=(17, 6))
# plot the x values vs the y values
ax[0].plot(x, y, color='purple', alpha=0.7, linewidth=0.3)
ax[0].set_title('x-y phase plane')
# plot the x values vs the z values
ax[1].plot(x, z, color='purple', alpha=0.7, linewidth=0.3)
ax[1].set_title('x-z phase plane')
# plot the y values vs the z values
ax[2].plot(y, z, color='purple', alpha=0.7, linewidth=0.3)
ax[2].set_title('y-z phase plane')
fig.savefig('{}/lorenz-attractor-phase-plane.png'.format(save_folder),
dpi=180, bbox_inches='tight')
plt.show()
The Manim code is incredibly extensive and it is perhaps a terrible idea to include it in this article — I will dedicate a separate one to it and explain the process. Meanwhile, find below a temporary render. I will do a more involved one as computation power permits but for now, this should suffice. It’s some pretty neat stuff! Shout out 3Blue1Brown for Manim.
That’s all from me here! I hope you enjoyed the article. As always, criticism is welcome.
Thank you for reading and have a great day!
