This is just a showcase of some cool matplotlib features, and also a visual treat 😎
We will, in various fashions, visualize a system of three ordinary differential equations that define the so-called Lorenz attractor:
Our figures will be defined by various parameters:
Attributions: The Lorenz function code is derived from the Matplotlib project website. The temporal animation code is entirely based on work by Ruben Esteche.
(Click blue links to navigate)
Part 1: Nested Plots
Part 2: Spatial Animation
Part 3: Temporal Animation
Appendix A: Offline code for spatial animation
Appendix B: Offline code for temporal animation
"From the beginning of chaos research until today, the unpredictability of chaos has been a central theme. It is widely believed and claimed by philosophers, mathematicians and physicists alike that chaos has a new implication for unpredictability, meaning that chaotic systems are unpredictable in a way that other deterministic systems are not. [...] However, this is not the case. [...] I will propose that the sought-after new implication of chaos for unpredictability is the following: for predicting any event, all sufficiently past events are approximately probabilistically irrelevant. [...] Chaotic behaviour is multi-faceted and takes various forms. Yet if the aim is to identify a general new implication of chaos for unpredictability, I think this is the best we can get."
Charlotte Werndl: What Are the New Implications of Chaos for Unpredictability? Brit. J. Phil. Sci. 60 (2009), 195–220)
Displayed below are 12 sets of respectively 3 solutions each.
Parameters are randomized for:
For non-mathematicians, this may facilitate developing an intuition for the interplay of the equations' parameters, and for the fundamental insight into the proximity of order and chaos.
# Import libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
display(HTML("<style>div.output_scroll { height: 999 }</style>"))
# Figure properties
fig = plt.figure(figsize=(10,120))
print("Lorenz Attractors (sigma, rho, beta )")
print("[ 3 partial solutions per subplot ]")
# Iterators and counters
plot_index = 0
for outer in range(12):
plot_index = plot_index + 1
# Projection parameters
is_pitch = 20
is_roll = 0
is_yaw = np.random.uniform(0, 90)
# Subplot properties
ax = fig.add_subplot(12, 1, plot_index, projection='3d')
ax.view_init(elev=is_pitch, azim=is_yaw, roll=is_roll)
ax.set(xticklabels=[])
ax.set(yticklabels=[])
ax.set(zticklabels=[])
ax.set_xlabel("X Axis", fontsize=7)
ax.set_ylabel("Y Axis", fontsize=7)
s_param = np.empty(3)
r_param = np.empty(3)
b_param = np.empty(3)
subplot_label = ("[" + (str(outer)) + "]: " + "Azimuth " + (str(int(is_yaw)) +"° ("))
color_list = ['blue', 'red', 'green']
series_label = ["", "", ""]
for inner in range (3):
s_param[inner]= np.random.uniform(5, 25)
r_param[inner] = np.random.uniform(5, 25)
b_param[inner] = np.random.uniform(0.1, 5)
# Lorenz function
def lorenz(xyz, *, s=s_param[inner], r=r_param[inner], b=b_param[inner]):
x, y, z = xyz
x_dot = s*(y - x)
y_dot = r*x - y - x*z
z_dot = x*y - b*z
return np.array([x_dot, y_dot, z_dot])
dt = 0.01
num_steps = 3000
xyzs = np.empty((num_steps + 1, 3)) # Need one more for the initial values
xyzs[0] = (0., 1., 1.05) # Set initial values
# Step through "time", calculating the partial derivatives at the current point
# and using them to estimate the next point
for i in range(num_steps):
xyzs[i + 1] = xyzs[i] + lorenz(xyzs[i]) * dt
# Plot partial solution
ax.plot(*xyzs.T, lw=0.5, color=color_list[inner])
# Append parameters for labeling
series_label[inner] = (color_list[inner] + ": " + str(int(s_param[inner])) + ", " +
str(int(r_param[inner])) + ", " + str(round(b_param[inner],1)) + " | ")
subplot_label = subplot_label + series_label[inner]
subplot_label = subplot_label[:-2] + ")"
ax.set_title(label=subplot_label, fontsize=10, color='DarkGrey')
#ax.plot(*xyzs.T, lw=0.5, color=(np.random.uniform(0, 0.7), np.random.uniform(0, 0.7), np.random.uniform(0, 0.7)))
plt.show()
Lorenz Attractors (sigma, rho, beta ) [ 3 partial solutions per subplot ]
A single solution for the parameters σ = 6 , ρ = 23 , β = 1 is animated in space, while the function plot is static in time.
The animation is not rendered inline, but displayed as a movie, prepared offline. The Python code for the animated plot can be found in Appendix A.
from IPython.display import HTML
HTML("""
<div align="middle">
<video width="90%" controls autoplay loop muted>
<source src="videos/Lorenz_animation.mov" type="video/mp4">
</video></div>""")
A single solution for the classic parameters σ = 10 , ρ = 28 , β = 8/3 (that Lorenz himself used in 1963) is animated in time, while the function plot is static in space.
The animation is not rendered inline, but displayed as an animated gif, prepared offline. The Python code for the constituent static images can be found in Appendix B.
import IPython.display as IPdisplay
IPdisplay.Image(url='videos/Lorenz_animate.gif')
# Imports
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np
# Define figure
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(projection='3d')
# Lorenz Attractor
def lorenz(xyz, *, s=6, r=23, b=1):
"""
Parameters '.
----------
xyz : array-like, shape (3,)
Point of interest in three-dimensional space.
s, r, b : float
Parameters defining the Lorenz attractor.
Returns
-------
xyz_dot : array, shape (3,)
Values of the Lorenz attractor's partial derivatives at *xyz*.
"""
x, y, z = xyz
x_dot = s*(y - x)
y_dot = r*x - y - x*z
z_dot = x*y - b*z
return np.array([x_dot, y_dot, z_dot])
dt = 0.01
num_steps = 10000
xyzs = np.empty((num_steps + 1, 3)) # Need one more for the initial values
xyzs[0] = (0., 1., 1.05) # Set initial values
# Step through "time", calculating the partial derivatives at the current point
# and using them to estimate the next point
for i in range(num_steps):
xyzs[i + 1] = xyzs[i] + lorenz(xyzs[i]) * dt
ax.plot(*xyzs.T, lw=0.75, color='b')
# Set the axis labels
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# This will hide the grids and backplanes
#ax.set_axis_off()
# Rotate the axes and update
for angle in range(0, 360*4 + 1):
# Normalize the angle to the range [-180, 180] for display
angle_norm = (angle + 180) % 360 - 180
# Cycle through a full rotation of elevation, then azimuth, roll, and all
elev = azim = roll = 0
if angle <= 360:
azim = angle_norm
#elif angle <= 360*2:
#elev = angle_norm
#elif angle <= 360*3:
#roll = angle_norm
else:
elev = azim = roll = angle_norm
# Update the axis view and title
ax.view_init(elev, azim, roll)
plt.title('Elevation: %d°, Azimuth: %d°, Roll: %d°' % (elev, azim, roll))
plt.draw()
# .02 = meditative; .001 = limit of performance
plt.pause(.02)
This link will bring you back to the Table of contents.
### Imports
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt, glob, os
import IPython.display as IPdisplay, matplotlib.font_manager as fm
from scipy.integrate import odeint
from mpl_toolkits.mplot3d.axes3d import Axes3D
from PIL import Image
### Saving to file
save_folder = 'images/lorenz-animate'
if not os.path.exists(save_folder):
os.makedirs(save_folder)
### Define the initial system state, initial parameters
### and the time points to solve for,
### evenly spaced between the start and end times
# (x, y, z positions in space)
initial_state = [0.1, 0, 0]
# sigma, rho, and beta
sigma = 10.
rho = 28.
beta = 8./3.
# time points
start_time = 1
end_time = 60
interval = 100
time_points = np.linspace(start_time, end_time, end_time * interval)
### Now for the Lorenz system:
def lorenz_system(current_state, t):
x, y, z = current_state
dx_dt = sigma * (y - x)
dy_dt = x * (rho - z) - y
dz_dt = x * y - beta * z
return [dx_dt, dy_dt, dz_dt]
### Plot the system line simplified in three dimensions
def plot_lorenz(xyz, n):
fig = plt.figure(figsize=(16, 12))
#############################################
# ORIGINAL CODE WAS USING DEPRECATED SYNTAX #
# ax = fig.gca(projection='3d')
# THIS IS THE WORKING SOLUTION BY m.e.: #
ax = fig.add_subplot(projection='3d')
#############################################
ax.xaxis.set_pane_color((1,1,1,1))
ax.yaxis.set_pane_color((1,1,1,1))
ax.zaxis.set_pane_color((1,1,1,1))
x = xyz[:, 0]
y = xyz[:, 1]
z = xyz[:, 2]
ax.plot(x, y, z, color='b', alpha=0.7, linewidth=0.7)
ax.set_xlim((-30,30))
ax.set_ylim((-30,30))
ax.set_zlim((0,50))
ax.set_title('Lorenz system (sigma=10, rho=28, beta=8/3)')
plt.savefig('{}/{:03d}.png'.format(save_folder, n), dpi=300, bbox_inches='tight', pad_inches=0.1)
plt.close()
### return a list in iteratively larger chunks,
### and handle them to facilitate future animation
def get_chunks(full_list, size):
size = max(1, size)
chunks = [full_list[0:i] for i in range(1, len(full_list) + 1, size)]
return chunks
chunks = get_chunks(time_points, size=20)
# get the points to plot, one chunk of time steps at a time,
# by integrating the system of equations
points = [odeint(lorenz_system, initial_state, chunk) for chunk in chunks]
# plot each set of points, one at a time, saving each plot
for n, point in enumerate(points):
plot_lorenz(point, n)
#######################################
# Somehow, GIF-creation did NOT work. #
# Thus I used a third party tool #
# to create a gif from the 300 PNGs, #
# which were fine! #
#######################################
'''
### Create an animated gif of all the plots
### then display it inline
# create a tuple of display durations, one for each frame
first_last = 100 #show the first and last frames for 100 ms
standard_duration = 5 #show all other frames for 5 ms
durations = tuple([first_last] + [standard_duration] * (len(points) - 2) + [first_last])
# load all the static images into a list
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:])
# verify that the number of frames in the gif equals the number of image files and durations
Image.open(gif_filepath).n_frames == len(images) == len(durations)
IPdisplay.Image(url=videos/Lorenz_animated.gif)
'''
#######################################
print("Done with preparation of PNGs!")
This link will bring you back to the Table of contents.