6.2 Animations

How we do visualization

Pre-computation Strategy

For computationally intensive tasks, we separate computation from visualization by pre-computing all required data and storing it in arrays. This strategy enables efficient animation through matplotlib’s FuncAnimation class, which simply loops through the pre-computed data to update the display. The same pre-computed arrays can also drive interactive exploration through sliders in marimo or Jupyter notebook widgets, allowing users to navigate through the solution space without recalculating at each step. This separation proves essential for problems like mechanism kinematics or chaotic attractors, where solving constraint equations or integrating differential equations at every frame would make smooth animation impossible.

Example: Three-Bar Linkage

This example shows how to animate a three-bar linkage mechanism using matplotlib’s FuncAnimation. The mechanism consists of three rigid bars connected by revolute joints, where we solve the kinematic constraints numerically to find valid configurations.

Code

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 60

# Mechanism parameters
r_OA = 100.0  # Length of link OA
r_CB = 75.0   # Length of link CB
rr_OC = np.array([250, 50])  
r_AB = np.sqrt((250 - 75)**2 + (100 - 50)**2)

def equations(vars, theta):
    """
    Constraint equations for three-bar linkage.
    
    Args:
        vars: [phi, gamma] - angles of links OA and AB in radians
        theta: Angle of link CB measured from vertical
    
    Returns:
        Residual vector that should be zero for valid configuration
    """
    phi, gamma = vars

    rr_OA = r_OA * np.array([np.sin(phi), np.cos(phi)]) 
    rr_AB = r_AB * np.array([np.sin(gamma), np.cos(gamma)])
    rr_OAB = rr_OA + rr_AB  # Position of B from chain O->A->B

    rr_CB = r_CB * np.array([-np.sin(theta), np.cos(theta)])
    rr_OCB = rr_OC + rr_CB  # Position of B from chain C->B

    return rr_OAB - rr_OCB

def is_valid_solution(theta, initial_guess, tolerance=1e-3):
    """Check if theta has a valid kinematic solution"""
    try:
        solution, info, ier, msg = scipy.optimize.fsolve(
            equations, initial_guess, args=(theta,), full_output=True
        )
        residual = np.linalg.norm(info['fvec'])
        return ier == 1 and residual < tolerance, solution
    except:
        return False, None

# Find valid theta range
print("Finding valid theta range for three-bar linkage...")
theta_test = np.linspace(np.deg2rad(-30), np.deg2rad(300), 330)
thetas = []
initial_guess = np.deg2rad([50, 70])

for theta in theta_test:
    is_valid, sol = is_valid_solution(theta, initial_guess)
    if is_valid:
        thetas.append(theta)
        initial_guess = sol

print(f"Valid theta range: [{np.rad2deg(thetas[0]):.1f}°, {np.rad2deg(thetas[-1]):.1f}°]")
print(f"Number of configurations: {len(thetas)}")

def plot_mechanism(ax, rr_OA, rr_OB, rr_OC, theta):
    """Plot the three-bar linkage"""
    ax.plot(*zip([0,0], rr_OA), '-ko', lw=2, label='Link OA')
    ax.plot(*zip(rr_OC, rr_OB), '-ro', lw=2, label='Link CB')
    ax.plot(*zip(rr_OA, rr_OB), '-bo', lw=2, label='Link AB')

    ax.text(rr_OA[0] + 5, rr_OA[1], 'A', fontsize=12, fontweight='bold')
    ax.text(rr_OB[0] + 5, rr_OB[1], 'B', fontsize=12, fontweight='bold')
    ax.text(rr_OC[0] + 5, rr_OC[1], 'C', fontsize=12, fontweight='bold')
    ax.text(5, 5, 'O', fontsize=12, fontweight='bold')
    
    ax.set_xlabel('x [mm]')
    ax.set_ylabel('y [mm]')
    ax.set_title(f'Three-Bar Linkage (θ = {np.rad2deg(theta):.1f}°)')
    ax.grid(True, alpha=0.3)
    ax.set_aspect('equal', adjustable='box')
    ax.set_xlim(-50, 300)
    ax.set_ylim(-50, 150)
    ax.legend()

fig, ax = plt.subplots(figsize=(10, 8))

def animate(frame_idx):
    """Animation function for each frame"""
    ax.clear()
    
    theta = thetas[frame_idx]
    initial_guess = np.deg2rad([50, 70]) 
    solution_rad = scipy.optimize.fsolve(equations, initial_guess, args=(theta,))
    phi_sol, gamma_sol = solution_rad

    rr_OA = r_OA * np.array([np.sin(phi_sol), np.cos(phi_sol)])
    rr_OB = rr_OA + r_AB * np.array([np.sin(gamma_sol), np.cos(gamma_sol)])

    plot_mechanism(ax, rr_OA, rr_OB, rr_OC, theta)

plt.close(fig)

# Create animation
anim = FuncAnimation(fig, animate, frames=len(thetas), 
                     interval=50, repeat=True)

# Convert to HTML5 video for display
video_html = anim.to_html5_video()
video_html = video_html.replace('controls', 'controls autoplay loop muted')
video_html = video_html.replace('<video ', '<video style="width: 100%; height: auto;" ')
HTML(video_html)

Finding valid theta range for three-bar linkage...
Valid theta range: [-1.9°, 204.7°]
Number of configurations: 207

Example: Crank-Slider Mechanism

The crank-slider mechanism converts rotational motion to linear motion (or vice versa). This is commonly found in piston engines and many mechanical systems. We solve the kinematic constraints to animate the mechanism.

Code

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 60

import warnings
warnings.filterwarnings('ignore', message='The iteration is not making good progress')

# Mechanism parameters
r = 50.0   # Crank length
L = 150.0  # Connecting rod length

def constraint_equations(vars, theta):
    """
    Constraint equations for crank-slider mechanism.
    
    Args:
        vars: [phi, x_A] - angle of connecting rod and slider position
        theta: crank angle (input)
        
    Returns:
        residuals: constraint equation residuals (should be zero)
    """
    phi, x_A = vars
    
    # Position of B (end of crank)
    rr_OB = r * np.array([-np.cos(theta), np.sin(theta)])
    
    # Position of A from B using connecting rod
    rr_BA = L * np.array([np.cos(phi), np.sin(phi)])
    
    # Position of A (computed from chain O->B->A)
    rr_OA_computed = rr_OB + rr_BA
    
    # Position of A (constrained to horizontal axis)
    rr_OA_constraint = np.array([x_A, 0])
    
    # Residuals (should be zero when constraints are satisfied)
    return rr_OA_computed - rr_OA_constraint

# Prepare animation data
theta_range = np.linspace(0, 2*np.pi, 60)
solutions = []
initial_guess = [np.deg2rad(180), -150]

print("Computing crank-slider configurations...")
for theta in theta_range:
    sol = scipy.optimize.fsolve(constraint_equations, initial_guess, args=(theta,))
    solutions.append((theta, sol[0], sol[1]))  # (theta, phi, x_A)
    initial_guess = sol

print(f"Computed {len(solutions)} configurations")

fig, ax = plt.subplots(figsize=(12, 6))

def plot_mechanism(theta, phi, x_A):
    """Plot the crank-slider mechanism for given configuration."""
    ax.clear()
    
    # Calculate positions
    O = np.array([0, 0])
    B = r * np.array([-np.cos(theta), np.sin(theta)])
    A = np.array([x_A, 0])
    
    # Plot ground
    ground_x = np.linspace(-250, 50, 100)
    ax.plot(ground_x, np.zeros_like(ground_x), 'k-', linewidth=3, alpha=0.3)
    
    # Plot slider rail
    ax.plot([-250, 50], [0, 0], 'k-', linewidth=5)
    
    # Plot crank OB
    ax.plot([O[0], B[0]], [O[1], B[1]], 'b-o', linewidth=3, markersize=8, label='Crank')
    
    # Plot connecting rod BA
    ax.plot([B[0], A[0]], [B[1], A[1]], 'r-o', linewidth=3, markersize=8, label='Connecting rod')
    
    # Plot slider
    slider_width = 20
    slider_height = 15
    slider_rect = plt.Rectangle((x_A - slider_width/2, -slider_height/2), 
                                slider_width, slider_height, 
                                facecolor='green', edgecolor='black', linewidth=2, label='Slider')
    ax.add_patch(slider_rect)
    
    # Plot origin
    ax.plot(O[0], O[1], 'ko', markersize=10)
    
    # Labels
    ax.text(B[0] + 5, B[1] + 5, 'B', fontsize=12, fontweight='bold')
    ax.text(A[0], A[1] + 10, 'A', fontsize=12, fontweight='bold')
    ax.text(O[0] - 10, O[1] + 10, 'O', fontsize=12, fontweight='bold')
    
    # Settings
    ax.set_xlim(-220, 60)
    ax.set_ylim(-80, 80)
    ax.set_aspect('equal', adjustable='box')
    ax.grid(True, alpha=0.3)
    ax.set_xlabel('x [mm]')
    ax.set_ylabel('y [mm]')
    ax.set_title(f'Crank-Slider Mechanism (θ = {np.rad2deg(theta):.1f}°, φ = {np.rad2deg(phi):.1f}°)')
    ax.legend(loc='upper right')

def animate(frame_idx):
    """Animation function."""
    theta, phi, x_A = solutions[frame_idx]
    plot_mechanism(theta, phi, x_A)

# Create animation
plt.close(fig)
anim = FuncAnimation(fig, animate, frames=len(solutions), 
                     interval=50, repeat=True)

# Convert to HTML5 video
video_html = anim.to_html5_video()
video_html = video_html.replace('controls', 'controls autoplay loop muted')
video_html = video_html.replace('<video ', '<video style="width: 100%; height: auto;" ')
HTML(video_html)

Computing crank-slider configurations...
Computed 60 configurations

Example: Simple Wave

A basic example showing how to animate a traveling sine wave. This demonstrates the core concept of matplotlib animations without complex physics.

Code

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 60

# Create x values
x = np.linspace(0, 2*np.pi, 200)

# Create figure and axis
fig, ax = plt.subplots(figsize=(10, 4))
line, = ax.plot([], [], 'b-', linewidth=2)

# Set up the plot limits
ax.set_xlim(0, 2*np.pi)
ax.set_ylim(-1.5, 1.5)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Traveling Sine Wave')
ax.grid(True, alpha=0.3)

def init():
    """Initialize animation"""
    line.set_data([], [])
    return line,

def animate(frame):
    """Animation function - updates for each frame"""
    # Calculate phase shift based on frame number
    phase = frame * 0.1
    
    # Calculate y values for traveling wave
    y = np.sin(x - phase)
    
    # Update the line data
    line.set_data(x, y)
    
    # Update title with current phase
    ax.set_title(f'Traveling Sine Wave (phase = {phase:.2f})')
    
    return line,

# Create animation
anim = FuncAnimation(fig, animate, init_func=init, 
                     frames=100, interval=50, blit=True, repeat=True)

plt.close(fig)

# Convert to HTML5 video
video_html = anim.to_html5_video()
video_html = video_html.replace('controls', 'controls autoplay loop muted')
video_html = video_html.replace('<video ', '<video style="width: 100%; height: auto;" ')
HTML(video_html)

Example: Lissajous Curve

Lissajous curves are parametric plots that create beautiful patterns. This example shows how to animate the drawing of a Lissajous curve by gradually revealing more of the curve over time.

Code

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 60

# Lissajous curve parameters
A = 1.0      # Amplitude in x
B = 1.0      # Amplitude in y
a = 3        # Frequency in x
b = 2        # Frequency in y
delta = np.pi/2  # Phase shift

# Create parameter t
t_full = np.linspace(0, 2*np.pi, 500)
x_full = A * np.sin(a * t_full + delta)
y_full = B * np.sin(b * t_full)

# Create figure
fig, ax = plt.subplots(figsize=(8, 8))
line, = ax.plot([], [], 'b-', linewidth=2, alpha=0.6)
point, = ax.plot([], [], 'ro', markersize=10)

ax.set_xlim(-1.5, 1.5)
ax.set_ylim(-1.5, 1.5)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title(f'Lissajous Curve (a={a}, b={b}, δ={delta:.2f})')
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')

def init():
    """Initialize animation"""
    line.set_data([], [])
    point.set_data([], [])
    return line, point

def animate(frame):
    """Animation function - gradually reveals the curve"""
    # Number of points to show
    n_points = int((frame / 100) * len(t_full))
    
    if n_points > 0:
        # Update the line (trail)
        line.set_data(x_full[:n_points], y_full[:n_points])
        
        # Update the point (current position)
        point.set_data([x_full[n_points-1]], [y_full[n_points-1]])
    
    return line, point

# Create animation
anim = FuncAnimation(fig, animate, init_func=init,
                     frames=100, interval=50, blit=True, repeat=True)

plt.close(fig)

# Convert to HTML5 video
video_html = anim.to_html5_video()
video_html = video_html.replace('controls', 'controls autoplay loop muted')
video_html = video_html.replace('<video ', '<video style="width: 100%; height: auto;" ')
HTML(video_html)

Example: Real-Time Interactive with PyQtGraph

For real-time interactive animations requiring user controls like sliders, buttons, and mouse interaction, PyQtGraph provides significantly better performance than matplotlib’s FuncAnimation. The Pendulum Simulation notebook demonstrates this approach, combining PyQt6 GUI components with PyQtGraph’s fast plotting capabilities. A QTimer drives the animation loop while signal/slot connections handle interactive events such as slider adjustments to modify system parameters, mouse clicks to set initial conditions, and keyboard controls for pause and reset functionality. Multiple synchronized plots update simultaneously, showing both the pendulum’s motion and its energy evolution over time. This architecture suits simulations requiring high frame rates above 60 fps with real-time parameter adjustment during execution, making it ideal for interactive physics demonstrations where users explore system behavior by manipulating parameters and observing immediate responses.

Example: Particle System

An animated scatter plot showing multiple particles moving and bouncing within a box. This demonstrates how to animate scatter plots with changing positions and colors.

Code

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 60

# Simulation parameters
n_particles = 50
box_size = 10.0
dt = 0.1

# Initialize particle positions and velocities
np.random.seed(42)
positions = np.random.uniform(0, box_size, (n_particles, 2))
velocities = np.random.uniform(-1, 1, (n_particles, 2))

# Create colors based on speed
def get_colors(velocities):
    """Calculate colors based on particle speed"""
    speeds = np.linalg.norm(velocities, axis=1)
    return speeds

# Create figure
fig, ax = plt.subplots(figsize=(8, 8))
scatter = ax.scatter([], [], s=100, c=[], cmap='viridis', vmin=0, vmax=2)
colorbar = plt.colorbar(scatter, ax=ax, label='Speed')

ax.set_xlim(0, box_size)
ax.set_ylim(0, box_size)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Particle System with Bouncing')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)

def init():
    """Initialize animation"""
    scatter.set_offsets(np.empty((0, 2)))
    return scatter,

def update_particles(positions, velocities, dt):
    """Update particle positions and handle wall bounces"""
    # Update positions
    positions += velocities * dt
    
    # Bounce off walls (x-direction)
    mask_x_low = positions[:, 0] < 0
    mask_x_high = positions[:, 0] > box_size
    velocities[mask_x_low, 0] *= -1
    velocities[mask_x_high, 0] *= -1
    positions[mask_x_low, 0] = -positions[mask_x_low, 0]
    positions[mask_x_high, 0] = 2*box_size - positions[mask_x_high, 0]
    
    # Bounce off walls (y-direction)
    mask_y_low = positions[:, 1] < 0
    mask_y_high = positions[:, 1] > box_size
    velocities[mask_y_low, 1] *= -1
    velocities[mask_y_high, 1] *= -1
    positions[mask_y_low, 1] = -positions[mask_y_low, 1]
    positions[mask_y_high, 1] = 2*box_size - positions[mask_y_high, 1]
    
    return positions, velocities

def animate(frame):
    """Animation function"""
    global positions, velocities
    
    # Update physics
    positions, velocities = update_particles(positions, velocities, dt)
    
    # Update scatter plot
    scatter.set_offsets(positions)
    scatter.set_array(get_colors(velocities))
    
    return scatter,

# Create animation
anim = FuncAnimation(fig, animate, init_func=init,
                     frames=100, interval=50, blit=True, repeat=True)

plt.close(fig)

# Convert to HTML5 video
video_html = anim.to_html5_video()
video_html = video_html.replace('controls', 'controls autoplay loop muted')
video_html = video_html.replace('<video ', '<video style="width: 100%; height: auto;" ')
HTML(video_html)

Example: Newton Fractal

We examine the convergence behavior of the Newton-Raphson method when applied to the polynomial equation \(z^5 + z^2 - z + 1\) across the complex plane. The method iterates according to \(z_{n+1} = z_n - f(z_n)/f'(z_n)\), where \(f(z) = z^5 + z^2 - z + 1\). This polynomial has five roots.

When we apply Newton’s method starting from different points in the complex plane, each initial condition converges to one of these five roots. We visualize this by coloring each pixel according to which root it ultimately reaches. Points that converge to the same root form a basin of attraction for that root. The fractal emerges at the boundaries between these basins, where infinitesimal changes in starting position lead to convergence toward entirely different roots.

This sensitivity to initial conditions characterizes chaotic dynamical systems. The boundaries exhibit self-similar structure at all scales: zooming into any boundary region reveals the same intricate patterns repeating infinitely. The animation demonstrates this property by zooming exponentially from a wide view into fine detail, then returning to the original scale.

This visualization connects numerical analysis (the Newton-Raphson method), complex dynamics (basin boundaries and attractors), and fractal geometry (self-similar structure across scales). The pre-computation strategy separates the expensive fractal calculation from the animation rendering, enabling smooth playback despite the computational intensity of iterating Newton’s method across hundreds of thousands of complex points for each frame.

See Grant Sanderson’s (3Blue1Brown) video explanation for more details on the mathematics and dynamics behind Newton fractals.

Code

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 60


import dill
with open('Newton_Fractals_Frames.dill', 'rb') as f:
    (fractal_images, zoom_levels, n_frames, frames_zoom_in, frames_zoom_out, 
     zoom_start, zoom_end, center_x, center_y, roots, ROOT_COLORS, n_roots) = dill.load(f)

# Create figure
fig, ax = plt.subplots(figsize=(10, 10))
fig.patch.set_facecolor('#222222')
ax.set_facecolor('#222222')

# Create colormap from root colors
from matplotlib.colors import ListedColormap
cmap = ListedColormap(ROOT_COLORS)

# Initialize image
img = ax.imshow(fractal_images[0], extent=[center_x - zoom_start, center_x + zoom_start,
                                            center_y - zoom_start, center_y + zoom_start],
                cmap=cmap, origin='lower', interpolation='bilinear', vmin=0, vmax=n_roots-1)

# Plot roots as markers with white rings
root_scatter = ax.scatter(roots.real, roots.imag, s=200, c='white', 
                          edgecolors='white', linewidths=3, zorder=10, alpha=0.8)
root_centers = ax.scatter(roots.real, roots.imag, s=50, c=ROOT_COLORS, 
                          edgecolors='white', linewidths=1, zorder=11)

ax.set_xlim(center_x - zoom_start, center_x + zoom_start)
ax.set_ylim(center_y - zoom_start, center_y + zoom_start)
ax.set_aspect('equal')
ax.axis('off')
ax.set_title('Newton Fractal: Basins of Attraction for $z^5 = 1$', 
             color='white', fontsize=14, pad=20)

def animate(frame):
    """Animation function - zoom in then out"""
    zoom = zoom_levels[frame]
    
    # Update fractal image
    img.set_data(fractal_images[frame])
    img.set_extent([center_x - zoom, center_x + zoom,
                   center_y - zoom, center_y + zoom])
    
    # Update axis limits
    ax.set_xlim(center_x - zoom, center_x + zoom)
    ax.set_ylim(center_y - zoom, center_y + zoom)
    
    # Update title with zoom level
    if frame < frames_zoom_in:
        phase = "Zooming in"
    else:
        phase = "Zooming out"
    ax.set_title(f'Newton Fractal: $z^5 = 1$ - {phase} (zoom: {1/zoom:.1f}x)', 
                color='white', fontsize=14, pad=20)
    
    return [img]

# Create animation with smooth 16ms interval (60fps)
# 1200 frames * 16ms = 19.2 seconds total 
fps = 60
interval = int(1000 / fps)
anim = FuncAnimation(fig, animate, frames=len(zoom_levels),
                     interval=interval, blit=True, repeat=True)

plt.close(fig)

# Convert to HTML5 video
video_html = anim.to_html5_video()
video_html = video_html.replace('controls', 'controls autoplay loop muted')
video_html = video_html.replace('<video ', '<video style="width: 100%; height: auto;" ')
HTML(video_html)

Example: Lorenz Attractor (3D)

The Lorenz attractor is one of the most famous chaotic dynamical systems, originally derived from atmospheric convection equations. It creates the iconic “butterfly” shape in 3D space and demonstrates sensitive dependence on initial conditions.

Code

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 60
from scipy.integrate import odeint

# Lorenz attractor parameters
sigma = 10.0
rho = 28.0
beta = 8.0 / 3.0

def lorenz_system(state, t):
    """
    Lorenz attractor differential equations.
    
    dx/dt = σ(y - x)
    dy/dt = x(ρ - z) - y
    dz/dt = xy - βz
    """
    x, y, z = 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]

# Generate 10 different initial conditions
np.random.seed(42)
n_trajectories = 10
initial_conditions = []
for i in range(n_trajectories):
    # Small random perturbations around (0, 1, 1.05)
    ic = [np.random.uniform(-1, 1), 
          np.random.uniform(0, 2), 
          np.random.uniform(20, 25)]
    initial_conditions.append(ic)

# Time array
t = np.linspace(0, 50, 5000)
dt_per_point = t[1] - t[0]

# Solve ODE for all initial conditions
# print("Computing Lorenz attractor trajectories...")
all_trajectories = []
for i, ic in enumerate(initial_conditions):
    trajectory = odeint(lorenz_system, ic, t)
    all_trajectories.append(trajectory)
    # print(f"Computed trajectory {i+1}/{n_trajectories}")

# Create figure with 3D axis
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')

# Set dark background
fig.patch.set_facecolor('#222222')
ax.set_facecolor('#222222')

# Make panes transparent
ax.xaxis.pane.fill = False
ax.yaxis.pane.fill = False
ax.zaxis.pane.fill = False

# Make pane edges invisible
ax.xaxis.pane.set_edgecolor('none')
ax.yaxis.pane.set_edgecolor('none')
ax.zaxis.pane.set_edgecolor('none')

# Colors for different trajectories using viridis colormap
cmap = plt.get_cmap('viridis')
colors = [cmap(i / n_trajectories) for i in range(n_trajectories)]

# Initialize lines and points for each trajectory
lines = []
points = []
for i in range(n_trajectories):
    line, = ax.plot([], [], [], '-', linewidth=2, alpha=0.8, color=colors[i])
    point, = ax.plot([], [], [], 'o', markersize=3, color=colors[i])
    lines.append(line)
    points.append(point)

# Find global limits across all trajectories
all_x = np.concatenate([traj[:, 0] for traj in all_trajectories])
all_y = np.concatenate([traj[:, 1] for traj in all_trajectories])
all_z = np.concatenate([traj[:, 2] for traj in all_trajectories])

x_min, x_max = all_x.min(), all_x.max()
y_min, y_max = all_y.min(), all_y.max()
z_min, z_max = all_z.min(), all_z.max()

# Set symmetric limits around center
x_center = (x_max + x_min) / 2
y_center = (y_max + y_min) / 2
z_center = (z_max + z_min) / 2

x_range = max(abs(x_max - x_center), abs(x_min - x_center))
y_range = max(abs(y_max - y_center), abs(y_min - y_center))
z_range = max(abs(z_max - z_center), abs(z_min - z_center))
max_range = max(x_range, y_range, z_range) * 1.1

ax.set_xlim(x_center - max_range, x_center + max_range)
ax.set_ylim(y_center - max_range, y_center + max_range)
ax.set_zlim(z_center - max_range, z_center + max_range)

# Remove all ticks
ax.set_xticks([])
ax.set_yticks([])
ax.set_zticks([])

# Remove grid
ax.grid(False)

# Hide the default axes completely
ax.xaxis.line.set_color('none')
ax.yaxis.line.set_color('none')
ax.zaxis.line.set_color('none')

# Hide axis labels from default axes
ax.set_xlabel('')
ax.set_ylabel('')
ax.set_zlabel('')

# Draw custom axis lines through the center with arrows
axis_color = 'lightgray'
axis_alpha = 0.6
arrow_length = max_range * 0.95

# X-axis with arrow
ax.quiver(x_center, y_center, z_center, arrow_length, 0, 0, 
          color=axis_color, alpha=axis_alpha, arrow_length_ratio=0.15, 
          linewidth=1.5, zorder=0)
ax.quiver(x_center, y_center, z_center, -arrow_length, 0, 0, 
          color=axis_color, alpha=axis_alpha*0.5, arrow_length_ratio=0, 
          linewidth=1.5, zorder=0)

# Y-axis with arrow
ax.quiver(x_center, y_center, z_center, 0, arrow_length, 0, 
          color=axis_color, alpha=axis_alpha, arrow_length_ratio=0.15, 
          linewidth=1.5, zorder=0)
ax.quiver(x_center, y_center, z_center, 0, -arrow_length, 0, 
          color=axis_color, alpha=axis_alpha*0.5, arrow_length_ratio=0, 
          linewidth=1.5, zorder=0)

# Z-axis with arrow
ax.quiver(x_center, y_center, z_center, 0, 0, arrow_length, 
          color=axis_color, alpha=axis_alpha, arrow_length_ratio=0.15, 
          linewidth=1.5, zorder=0)
ax.quiver(x_center, y_center, z_center, 0, 0, -arrow_length, 
          color=axis_color, alpha=axis_alpha*0.5, arrow_length_ratio=0, 
          linewidth=1.5, zorder=0)

# Add axis labels at the end of each axis
label_offset = max_range * 1.05
ax.text(x_center + label_offset, y_center, z_center, 'X', color='lightgray', fontsize=12, ha='center')
ax.text(x_center, y_center + label_offset, z_center, 'Y', color='lightgray', fontsize=12, ha='center')
ax.text(x_center, y_center, z_center + label_offset, 'Z', color='lightgray', fontsize=12, ha='center')

ax.set_title('Lorenz Attractor', color='lightgray', fontsize=14, pad=20)

# Adjust viewing angle for better visualization
ax.view_init(elev=20, azim=45)

def init():
    """Initialize animation"""
    for line, point in zip(lines, points):
        line.set_data([], [])
        line.set_3d_properties([])
        point.set_data([], [])
        point.set_3d_properties([])
    return lines + points

def animate(frame):
    """Animation function - gradually reveals the trajectories with fading"""
    # Slower animation - show fewer points per frame
    n_points = int((frame / 200) * len(t))
    
    if n_points > 1:
        # Calculate current time
        current_time = t[min(n_points-1, len(t)-1)]
        
        # Calculate visible window based on fading rules
        # Start fading at t=10, complete fade by t=30
        if current_time < 10:
            # Before t=10: show everything from start
            start_idx = 0
        elif current_time <= 30:
            # From t=10 to t=30: gradually fade from beginning
            fade_progress = (current_time - 10) / 20  # 0 to 1
            start_time = fade_progress * (current_time - 10)
            start_idx = int(start_time / dt_per_point)
        else:
            # After t=30: maintain 20-second sliding window
            start_time = current_time - 20
            start_idx = int(start_time / dt_per_point)
        
        # Ensure indices are within bounds
        start_idx = max(0, start_idx)
        end_idx = min(n_points, len(t))
        
        # Update each trajectory
        for i, (line, point, traj) in enumerate(zip(lines, points, all_trajectories)):
            x_traj = traj[:, 0]
            y_traj = traj[:, 1]
            z_traj = traj[:, 2]
            
            # Update the trajectory line with visible window
            line.set_data(x_traj[start_idx:end_idx], y_traj[start_idx:end_idx])
            line.set_3d_properties(z_traj[start_idx:end_idx])
            
            # Update the current position point
            if end_idx > 0:
                point.set_data([x_traj[end_idx-1]], [y_traj[end_idx-1]])
                point.set_3d_properties([z_traj[end_idx-1]])
        
        # Rotate the view slightly for dynamic visualization
        ax.view_init(elev=20, azim=45 + frame * 0.3)
    
    return lines + points

# Create animation with more frames for slower motion
anim = FuncAnimation(fig, animate, init_func=init,
                     frames=200, interval=50, blit=False, repeat=True)

plt.close(fig)

# Convert to HTML5 video
video_html = anim.to_html5_video()
video_html = video_html.replace('controls', 'controls autoplay loop muted')
video_html = video_html.replace('<video ', '<video style="width: 100%; height: auto;" ')
HTML(video_html)

Example: Aizawa Attractor (3D)

The Aizawa attractor is a chaotic dynamical system that creates beautiful 3D trajectories. This example demonstrates 3D animation in matplotlib, showing how the system evolves over time with a trailing path.

Code

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 60
from scipy.integrate import odeint

# Aizawa attractor parameters
a = 0.95
b = 0.7
c = 0.6
d = 3.5
e = 0.25
f = 0.1

def aizawa_system(state, t):
    """
    Aizawa attractor differential equations.
    
    dx/dt = (z - b)*x - d*y
    dy/dt = d*x + (z - b)*y
    dz/dt = c + a*z - z³/3 - (x² + y²)
    """
    x, y, z = state
    
    dx_dt = (z - b) * x - d * y
    dy_dt = d * x + (z - b) * y
    dz_dt = c + a * z - (z**3) / 3 - (x**2 + y**2)
    
    return [dx_dt, dy_dt, dz_dt]

# Generate 10 different initial conditions
np.random.seed(42)
n_trajectories = 40
initial_conditions = []
for i in range(n_trajectories):
    # Small random perturbations around origin
    ic = [np.random.uniform(-0.2, 0.2), 
          np.random.uniform(-0.2, 0.2), 
          np.random.uniform(-0.2, 0.2)]
    initial_conditions.append(ic)

# Time array
t = np.linspace(0, 30, 5000)
dt_per_point = t[1] - t[0]

# Solve ODE for all initial conditions
# print("Computing Aizawa attractor trajectories...")
all_trajectories = []
for i, ic in enumerate(initial_conditions):
    trajectory = odeint(aizawa_system, ic, t)
    all_trajectories.append(trajectory)
    # print(f"Computed trajectory {i+1}/{n_trajectories}")

# Create figure with 3D axis
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')

# Set dark background
fig.patch.set_facecolor('#222222')
ax.set_facecolor('#222222')

# Make panes transparent
ax.xaxis.pane.fill = False
ax.yaxis.pane.fill = False
ax.zaxis.pane.fill = False

# Make pane edges invisible
ax.xaxis.pane.set_edgecolor('none')
ax.yaxis.pane.set_edgecolor('none')
ax.zaxis.pane.set_edgecolor('none')

# Colors for different trajectories using viridis colormap
cmap = plt.get_cmap('viridis')
colors = [cmap(i / n_trajectories) for i in range(n_trajectories)]

# Initialize lines and points for each trajectory
lines = []
points = []
for i in range(n_trajectories):
    line, = ax.plot([], [], [], '-', linewidth=2, alpha=0.8, color=colors[i])
    point, = ax.plot([], [], [], 'o', markersize=3, color=colors[i])
    lines.append(line)
    points.append(point)

# Find global limits across all trajectories
all_x = np.concatenate([traj[:, 0] for traj in all_trajectories])
all_y = np.concatenate([traj[:, 1] for traj in all_trajectories])
all_z = np.concatenate([traj[:, 2] for traj in all_trajectories])

x_min, x_max = all_x.min(), all_x.max()
y_min, y_max = all_y.min(), all_y.max()
z_min, z_max = all_z.min(), all_z.max()

# Set symmetric limits around center
x_range = max(abs(x_min), abs(x_max))
y_range = max(abs(y_min), abs(y_max))
z_range = max(abs(z_min), abs(z_max))
max_range = max(x_range, y_range, z_range) * 1.1

ax.set_xlim(-max_range, max_range)
ax.set_ylim(-max_range, max_range)
ax.set_zlim(-max_range, max_range)

# Remove all ticks
ax.set_xticks([])
ax.set_yticks([])
ax.set_zticks([])

# Remove grid
ax.grid(False)

# Hide the default axes completely
ax.xaxis.line.set_color('none')
ax.yaxis.line.set_color('none')
ax.zaxis.line.set_color('none')

# Hide axis labels from default axes
ax.set_xlabel('')
ax.set_ylabel('')
ax.set_zlabel('')

# Draw custom axis lines through the center (origin) with arrows
axis_color = 'lightgray'
axis_alpha = 0.6
arrow_length = max_range * 0.95

# X-axis with arrow
ax.quiver(0, 0, 0, arrow_length, 0, 0, 
          color=axis_color, alpha=axis_alpha, arrow_length_ratio=0.15, 
          linewidth=1.5, zorder=0)
ax.quiver(0, 0, 0, -arrow_length, 0, 0, 
          color=axis_color, alpha=axis_alpha*0.5, arrow_length_ratio=0, 
          linewidth=1.5, zorder=0)

# Y-axis with arrow
ax.quiver(0, 0, 0, 0, arrow_length, 0, 
          color=axis_color, alpha=axis_alpha, arrow_length_ratio=0.15, 
          linewidth=1.5, zorder=0)
ax.quiver(0, 0, 0, 0, -arrow_length, 0, 
          color=axis_color, alpha=axis_alpha*0.5, arrow_length_ratio=0, 
          linewidth=1.5, zorder=0)

# Z-axis with arrow
ax.quiver(0, 0, 0, 0, 0, arrow_length, 
          color=axis_color, alpha=axis_alpha, arrow_length_ratio=0.15, 
          linewidth=1.5, zorder=0)
ax.quiver(0, 0, 0, 0, 0, -arrow_length, 
          color=axis_color, alpha=axis_alpha*0.5, arrow_length_ratio=0, 
          linewidth=1.5, zorder=0)

# Add axis labels at the end of each axis
label_offset = max_range * 1.05
ax.text(label_offset, 0, 0, 'X', color='lightgray', fontsize=12, ha='center')
ax.text(0, label_offset, 0, 'Y', color='lightgray', fontsize=12, ha='center')
ax.text(0, 0, label_offset, 'Z', color='lightgray', fontsize=12, ha='center')

ax.set_title('Aizawa Attractor', color='lightgray', fontsize=14, pad=20)

# Adjust viewing angle for better visualization
ax.view_init(elev=20, azim=45)

def init():
    """Initialize animation"""
    for line, point in zip(lines, points):
        line.set_data([], [])
        line.set_3d_properties([])
        point.set_data([], [])
        point.set_3d_properties([])
    return lines + points

def animate(frame):
    """Animation function - gradually reveals the trajectories with fading"""
    # Slower animation - show fewer points per frame
    n_points = int((frame / 200) * len(t))
    
    if n_points > 1:
        # Calculate current time
        current_time = t[min(n_points-1, len(t)-1)]
        
        # Calculate visible window based on fading rules
        # Start fading at t=10, complete fade by t=30
        t_anim_start = 15
        if current_time < t_anim_start:
            # Before t=10: show everything from start
            start_idx = 0
        elif current_time <= 30:
            # From t=10 to t=30: gradually fade from beginning
            # Window start moves from 0 to (current_time - 20)
            fade_progress = (current_time - t_anim_start) / 20  # 0 to 1
            start_time = fade_progress * (current_time - t_anim_start)
            start_idx = int(start_time / dt_per_point)
        else:
            # After t=30: maintain 20-second sliding window
            start_time = current_time - 20
            start_idx = int(start_time / dt_per_point)
        
        # Ensure indices are within bounds
        start_idx = max(0, start_idx)
        end_idx = min(n_points, len(t))
        
        # Update each trajectory
        for i, (line, point, traj) in enumerate(zip(lines, points, all_trajectories)):
            x_traj = traj[:, 0]
            y_traj = traj[:, 1]
            z_traj = traj[:, 2]
            
            # Update the trajectory line with visible window
            line.set_data(x_traj[start_idx:end_idx], y_traj[start_idx:end_idx])
            line.set_3d_properties(z_traj[start_idx:end_idx])
            
            # Update the current position point
            if end_idx > 0:
                point.set_data([x_traj[end_idx-1]], [y_traj[end_idx-1]])
                point.set_3d_properties([z_traj[end_idx-1]])
        
        # Rotate the view slightly for dynamic visualization
        ax.view_init(elev=20, azim=45 + frame * 0.3)
    
    return lines + points

# Create animation with more frames for slower motion
anim = FuncAnimation(fig, animate, init_func=init,
                     frames=200, interval=50, blit=False, repeat=True)

plt.close(fig)

# Convert to HTML5 video
video_html = anim.to_html5_video()
video_html = video_html.replace('controls', 'controls autoplay loop muted')
video_html = video_html.replace('<video ', '<video style="width: 100%; height: auto;" ')
HTML(video_html)

A note on visualization of mathematical concepts

Visualization of Mathematical Concepts

Mathematical visualization transforms abstract concepts into perceivable phenomena. When we animate the Newton fractal’s zoom sequence or watch the Lorenz attractor trace its butterfly wings, we engage with mathematics through perception, revealing structure that remains hidden in equations alone. This direct visual engagement builds intuition complementing analytical understanding, showing us how systems behave in addition to what equations describe.

Grant Sanderson’s 3Blue1Brown demonstrated the pedagogical power of high-quality mathematical animation, showing that visualization serves not as substitute for formal reasoning but as complementary mode of understanding. His work catalyzed a thriving community including channels like Welch Labs, sudgylacmoe, Polyamathematics, and Two Swaps, establishing visualization as legitimate pedagogical methodology. The Newton fractal example draws from this tradition, using the 3Blue1Brown color scheme to reveal self-similar fractal boundaries between basins of attraction that static images cannot capture.

The animations throughout this notebook reflect this philosophy. We do not merely show static plots but animate evolution over time, revealing dynamics our visual system grasps more immediately than symbolic descriptions alone permit. The three-bar linkage’s constrained motion, the chaotic attractors’ fading trajectories, the fractal boundaries emerging at finer scales—these animations transform abstract mathematical objects into phenomena we can observe and explore. If in doubt, write some visualization code.

Marimo notebooks

Check out python.ju.se/marimo as well as mechanics.ju.se for additional notebooks and study the code used to create these animations.

Animation Techniques and Patterns

Marimo sliders

The basic gist is to pre-compute the data and use marimo’s slider functionality to visualize the data interactively.

FuncAnimation Basics

The FuncAnimation class is the core of matplotlib animations. It works by repeatedly calling a function to update the plot. Similar to the marimo slider approach, we can pre-compute data for each frame and then use FuncAnimation to display it.

Basic structure:

from matplotlib.animation import FuncAnimation

fig, ax = plt.subplots()

def animate(frame):
    # Update plot for this frame
    # frame is the frame number (0, 1, 2, ...)
    pass

anim = FuncAnimation(fig, animate, frames=100, interval=50)

Key parameters: - fig: The figure object to animate - animate: Function called for each frame - frames: Number of frames (or iterable of frame data) - interval: Delay between frames in milliseconds - blit=True: Only redraw changed parts (faster) - repeat=True: Loop the animation - init_func: Function to initialize the animation

Pre-computing vs Real-time Updates

Pre-compute approach (recommended for complex calculations):

# Compute all data first
data = []
for i in range(n_frames):
    # Expensive computation here
    result = compute_frame(i)
    data.append(result)

# Then animate the pre-computed data
def animate(frame):
    # Just update the plot with pre-computed data
    line.set_data(x_data[frame], y_data[frame])

Real-time approach (for simple calculations):

def animate(frame):
    # Compute and display in the same function
    y = np.sin(x + frame * 0.1)
    line.set_data(x, y)

Converting Animations to Video

For Jupyter notebooks, convert animations to HTML5 video:

# Convert to HTML5 video
video_html = anim.to_html5_video()
video_html = video_html.replace('controls', 'controls autoplay loop muted')
video_html = video_html.replace('<video ', '<video style="width: 100%; height: auto;" ')
HTML(video_html)

Or use JavaScript-based animation:

HTML(anim.to_jshtml())

Performance Tips

Use blit=True when possible - only redraws changed elements
Pre-compute expensive calculations before animation
Close the figure with plt.close(fig) to prevent duplicate display
Limit frame rate with the interval parameter
Use appropriate data structures - arrays are faster than lists

Saving Animations

Save animations to video files:

# Save as MP4 (requires ffmpeg)
anim.save('animation.mp4', writer='ffmpeg', fps=30)

# Save as GIF (requires imagemagick or pillow)
anim.save('animation.gif', writer='pillow', fps=30)

Common Patterns

Pattern 1: Line plot animation

line, = ax.plot([], [], 'b-')

def animate(frame):
    line.set_data(x[:frame], y[:frame])
    return line,

Pattern 2: Scatter plot animation

scatter = ax.scatter([], [])

def animate(frame):
    scatter.set_offsets(positions[frame])
    scatter.set_array(colors[frame])
    return scatter,

Pattern 3: Multiple elements

line1, = ax.plot([], [], 'b-')
line2, = ax.plot([], [], 'r-')

def animate(frame):
    line1.set_data(x1[frame], y1[frame])
    line2.set_data(x2[frame], y2[frame])
    return line1, line2