Strange Attractor¶

Visualizations of the Lorenz System¶

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:

equations

Our figures will be defined by various parameters:

  • Equation parameters σ, ρ, β
  • Number of iterations through 'time'
  • Observer vantage (pitch, roll and yaw of graph)
  • Line properties (color, weight, translucency)

Attributions: The Lorenz function code is derived from the Matplotlib project website. The temporal animation code is entirely based on work by Ruben Esteche.

Table of Contents:¶

(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)

Part 1: Nested Plots¶

Displayed below are 12 sets of respectively 3 solutions each.

Parameters are randomized for:

  • yaw (rotation),
  • equations parameters σ, ρ, β.

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.

In [64]:
# 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 ]

Part 2: Spatial Animation¶

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.

In [1]:
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>""")
Out[1]:

This link will bring you back to the Table of contents.


Part 3: Temporal Animation¶

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.

In [2]:
import IPython.display as IPdisplay
IPdisplay.Image(url='videos/Lorenz_animate.gif')
Out[2]:

This link will bring you back to the Table of contents.


Appendix A:¶

Offline code for spatial animation¶


# 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.


Appendix B:¶

Offline code for temporal animation¶


### 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.


In [ ]: