Can a difference of just 0.001 radians (0.057 degrees) completely shatter determinism? In this investigation, we derive the equations of motion for the chaotic double pendulum from first principles using Lagrangian Mechanics, solve the non-linear coupled differential equations using a 4th-Order Runge-Kutta integrator, and animate the resulting 4D phase space attractor using Manim.
Figure 1: Numerical simulation of two coupled double pendulums differing by only Δθ₀ = 0.001 rad.
1. The Coordinate Frame and Kinematics
Consider two rigid, massless rods of lengths $L_1$ and $L_2$ holding point masses $m_1$ and $m_2$ in a uniform gravitational field $g$. We define the generalized coordinates as the deflection angles $\theta_1$ and $\theta_2$ measured relative to the vertical downward axis.
The Cartesian coordinates for the first bob $m_1$ are:
$$x_1 = L_1 \sin\theta_1, \quad y_1 = -L_1 \cos\theta_1$$
The Cartesian coordinates for the second bob $m_2$ depend on the position of the first:
$$x_2 = L_1 \sin\theta_1 + L_2 \sin\theta_2$$
$$y_2 = -L_1 \cos\theta_1 – L_2 \cos\theta_2$$
Differentiating each with respect to time $t$ via the chain rule yields the linear velocity components:
$$\dot{x}_1 = L_1 \dot{\theta}_1 \cos\theta_1, \quad \dot{y}_1 = L_1 \dot{\theta}_1 \sin\theta_1$$
$$\dot{x}_2 = L_1 \dot{\theta}_1 \cos\theta_1 + L_2 \dot{\theta}_2 \cos\theta_2$$
$$\dot{y}_2 = L_1 \dot{\theta}_1 \sin\theta_1 + L_2 \dot{\theta}_2 \sin\theta_2$$
$$v_1^2 = \dot{x}_1^2 + \dot{y}_1^2 = L_1^2 \dot{\theta}_1^2 (\cos^2\theta_1 + \sin^2\theta_1) = L_1^2 \dot{\theta}_1^2$$
$$v_2^2 = \dot{x}_2^2 + \dot{y}_2^2 = L_1^2 \dot{\theta}_1^2 + L_2^2 \dot{\theta}_2^2 + 2 L_1 L_2 \dot{\theta}_1 \dot{\theta}_2 \cos(\theta_1 – \theta_2)$$
2. Lagrangian Formulation ($L = T – V$)
The total kinetic energy $T$ of the mechanical system is the sum of the kinetic energies of both point masses:
$$T = \frac{1}{2} m_1 v_1^2 + \frac{1}{2} m_2 v_2^2$$
$$T = \frac{1}{2}(m_1 + m_2) L_1^2 \dot{\theta}_1^2 + \frac{1}{2} m_2 L_2^2 \dot{\theta}_2^2 + m_2 L_1 L_2 \dot{\theta}_1 \dot{\theta}_2 \cos(\theta_1 – \theta_2)$$
Taking $y=0$ at the fixed suspension pivot, the total gravitational potential energy $V$ is:
$$V = m_1 g y_1 + m_2 g y_2 = -(m_1 + m_2) g L_1 \cos\theta_1 – m_2 g L_2 \cos\theta_2$$
The Lagrangian $L = T – V$ is therefore:
$$L = \frac{1}{2}(m_1 + m_2) L_1^2 \dot{\theta}_1^2 + \frac{1}{2} m_2 L_2^2 \dot{\theta}_2^2 + m_2 L_1 L_2 \dot{\theta}_1 \dot{\theta}_2 \cos(\theta_1 – \theta_2) + (m_1 + m_2) g L_1 \cos\theta_1 + m_2 g L_2 \cos\theta_2$$
3. Euler-Lagrange Equations of Motion
According to Hamilton’s Principle of Stationary Action, the physical trajectories satisfy the Euler-Lagrange equations for each generalized coordinate $\theta_i$:
$$\frac{d}{dt}\left( \frac{\partial L}{\partial \dot{\theta}_1} \right) – \frac{\partial L}{\partial \theta_1} = 0, \qquad \frac{d}{dt}\left( \frac{\partial L}{\partial \dot{\theta}_2} \right) – \frac{\partial L}{\partial \theta_2} = 0$$
Evaluating the partial derivatives for $\theta_1$:
$$\frac{\partial L}{\partial \dot{\theta}_1} = (m_1 + m_2) L_1^2 \dot{\theta}_1 + m_2 L_1 L_2 \dot{\theta}_2 \cos(\theta_1 – \theta_2)$$
$$\frac{\partial L}{\partial \theta_1} = -m_2 L_1 L_2 \dot{\theta}_1 \dot{\theta}_2 \sin(\theta_1 – \theta_2) – (m_1 + m_2) g L_1 \sin\theta_1$$
Taking the total time derivative $\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta}_1}\right)$ and substituting gives the first equation of motion:
$$(m_1 + m_2) L_1 \ddot{\theta}_1 + m_2 L_2 \ddot{\theta}_2 \cos(\theta_1 – \theta_2) + m_2 L_2 \dot{\theta}_2^2 \sin(\theta_1 – \theta_2) + (m_1 + m_2) g \sin\theta_1 = 0$$
Repeating the exact procedure for $\theta_2$ gives the second coupled equation:
$$m_2 L_2 \ddot{\theta}_2 + m_2 L_1 \ddot{\theta}_1 \cos(\theta_1 – \theta_2) – m_2 L_1 \dot{\theta}_1^2 \sin(\theta_1 – \theta_2) + m_2 g \sin\theta_2 = 0$$
4. Solving for Angular Accelerations ($\ddot{\theta}_1, \ddot{\theta}_2$)
Notice that the system is linear in the second derivatives $\ddot{\theta}_1$ and $\ddot{\theta}_2$. We formulate this as a $2\times 2$ matrix system:
$$\begin{bmatrix} (m_1 + m_2) L_1 & m_2 L_2 \cos\Delta \\ L_1 \cos\Delta & L_2 \end{bmatrix} \begin{bmatrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{bmatrix} = \begin{bmatrix} -m_2 L_2 \dot{\theta}_2^2 \sin\Delta – (m_1 + m_2) g \sin\theta_1 \\ L_1 \dot{\theta}_1^2 \sin\Delta – g \sin\theta_2 \end{bmatrix}$$
where $\Delta = \theta_1 – \theta_2$.
Applying Cramer’s rule or direct matrix inversion yields the exact equations for $\ddot{\theta}_1$ and $\ddot{\theta}_2$ used in our simulation:
$$\ddot{\theta}_1 = \frac{-g(2m_1 + m_2)\sin\theta_1 – m_2 g \sin(\theta_1 – 2\theta_2) – 2\sin\Delta m_2 (\dot{\theta}_2^2 L_2 + \dot{\theta}_1^2 L_1 \cos\Delta)}{L_1 (2m_1 + m_2 – m_2 \cos(2\Delta))}$$
$$\ddot{\theta}_2 = \frac{2\sin\Delta (\dot{\theta}_1^2 L_1 (m_1 + m_2) + g(m_1 + m_2)\cos\theta_1 + \dot{\theta}_2^2 L_2 m_2 \cos\Delta)}{L_2 (2m_1 + m_2 – m_2 \cos(2\Delta))}$$
5. Fourth-Order Runge-Kutta (RK4) Integration
Because these non-linear ODEs have no closed-form analytical solution, we define the state vector $\mathbf{y}(t) = [\theta_1, \omega_1, \theta_2, \omega_2]^T$ and integrate forward in time using the classical 4th-Order Runge-Kutta method:
$$\mathbf{k}_1 = f(t_n, \mathbf{y}_n)$$
$$\mathbf{k}_2 = f\left(t_n + \frac{\Delta t}{2}, \mathbf{y}_n + \frac{\Delta t}{2}\mathbf{k}_1\right)$$
$$\mathbf{k}_3 = f\left(t_n + \frac{\Delta t}{2}, \mathbf{y}_n + \frac{\Delta t}{2}\mathbf{k}_2\right)$$
$$\mathbf{k}_4 = f(t_n + \Delta t, \mathbf{y}_n + \Delta t \mathbf{k}_3)$$
$$\mathbf{y}_{n+1} = \mathbf{y}_n + \frac{\Delta t}{6}(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4)$$
6. The Lyapunov Exponent & Sensitivity to Initial Conditions
In our Manim animation, we release two pendulums simultaneously:
- System A (Cyan): $\theta_1(0) = 90.000^\circ$, $\theta_2(0) = 90.000^\circ$
- System B (Electric Pink): $\theta_1(0) = 90.057^\circ$ (differing by only $0.001\text{ rad}$), $\theta_2(0) = 90.000^\circ$
For the first 3 seconds, both trajectories overlap so perfectly that they appear to be a single system. However, the system’s phase space separation $\|\Delta \mathbf{Z}(t)\|$ grows exponentially:
$$\|\Delta \mathbf{Z}(t)\| \approx \|\Delta \mathbf{Z}(0)\| e^{\lambda t}$$
where $\lambda > 0$ is the Maximal Lyapunov Exponent. By $t = 6\text{ seconds}$, the pink and cyan bobs are in completely opposite regions of space, demonstrating the quintessential hallmark of deterministic chaos: complete microscopic unpredictability emerging from fully deterministic laws.
7. Full Manim Community Source Code
Below is the complete, self-contained Python source code using Manim Community Edition. It includes both the 16:9 landscape master composition and the 9:16 vertical layout for YouTube Shorts and Instagram Reels.
Python 3
"""
Double Pendulum & Chaos Theory: From First Principles to Chaos
Engine: Manim Community Edition
Aesthetic: 3Blue1Brown Neo-Obsidian Dark Mode with Neon Cyan & Magenta Attractors
"""
import numpy as np
from manim import *
# Aesthetic Color Palette
COLOR_BG = "#0B0C10"
COLOR_ROD = "#C5C6C7"
COLOR_PIVOT = "#66FCF1"
COLOR_CYAN = "#45A29E"
COLOR_NEON_A = "#00F2FE" # First pendulum glow (Cyan)
COLOR_NEON_B = "#FF007F" # Second pendulum glow (Electric Pink)
COLOR_TEXT = "#EAEAEA"
COLOR_ACCENT = "#FFD166"
# Physics Parameters
G = 9.81 # Acceleration due to gravity (m/s^2)
L1 = 2.0 # Length of rod 1 (m)
L2 = 1.8 # Length of rod 2 (m)
M1 = 1.5 # Mass 1 (kg)
M2 = 1.0 # Mass 2 (kg)
def double_pendulum_derivatives(state, t):
"""
Computes [theta1_dot, omega1_dot, theta2_dot, omega2_dot] from Euler-Lagrange equations.
State vector: [th1, w1, th2, w2]
"""
th1, w1, th2, w2 = state
delta = th1 - th2
# Denominator terms
den1 = L1 * (2 * M1 + M2 - M2 * np.cos(2 * th1 - 2 * th2))
den2 = L2 * (2 * M1 + M2 - M2 * np.cos(2 * th1 - 2 * th2))
# Angular acceleration 1
num1 = (-G * (2 * M1 + M2) * np.sin(th1)
- M2 * G * np.sin(th1 - 2 * th2)
- 2 * np.sin(delta) * M2 * (w2**2 * L2 + w1**2 * L1 * np.cos(delta)))
alpha1 = num1 / den1
# Angular acceleration 2
num2 = (2 * np.sin(delta) * (w1**2 * L1 * (M1 + M2)
+ G * (M1 + M2) * np.cos(th1)
+ w2**2 * L2 * M2 * np.cos(delta)))
alpha2 = num2 / den2
return np.array([w1, alpha1, w2, alpha2])
def rk4_step(state, dt):
"""Fourth-order Runge-Kutta integrator for high numerical stability."""
k1 = double_pendulum_derivatives(state, 0)
k2 = double_pendulum_derivatives(state + 0.5 * dt * k1, 0)
k3 = double_pendulum_derivatives(state + 0.5 * dt * k2, 0)
k4 = double_pendulum_derivatives(state + dt * k3, 0)
return state + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
class DoublePendulumChaos(Scene):
"""
Landscape presentation (16:9) demonstrating sensitive dependence
on initial conditions (Butterfly Effect & Lyapunov divergence).
"""
def construct(self):
self.camera.background_color = COLOR_BG
# Title Card
title = Text("The Double Pendulum: Order to Chaos", font="Sans", font_size=38, color=COLOR_TEXT)
subtitle = Text("Euler-Lagrange Mechanics & The Butterfly Effect", font="Sans", font_size=22, color=COLOR_PIVOT)
VGroup(title, subtitle).arrange(DOWN, buff=0.2).to_edge(UP, buff=0.4)
self.play(FadeIn(title, shift=DOWN*0.3), FadeIn(subtitle, shift=UP*0.3), run_time=1.5)
self.wait(0.5)
# Scale and coordinate origin
SCALE = 0.95
pivot_point = UP * 1.5 + LEFT * 0.0
# Pivot anchor
pivot_dot = Dot(pivot_point, radius=0.12, color=COLOR_PIVOT)
pivot_ring = Circle(radius=0.22, color=COLOR_PIVOT, stroke_width=2).move_to(pivot_point)
self.play(Create(pivot_dot), Create(pivot_ring), run_time=1.0)
# Initial Conditions: Two systems differing by only 0.001 radians (0.057 degrees!)
# State: [theta1, omega1, theta2, omega2]
init_th1 = np.pi / 2.0 # 90 degrees horizontal
init_th2 = np.pi / 2.0
delta_init = 0.001 # 0.001 rad perturbation
state_A = np.array([init_th1, 0.0, init_th2, 0.0], dtype=float)
state_B = np.array([init_th1 + delta_init, 0.0, init_th2, 0.0], dtype=float)
# Setup Visual Elements for System A (Cyan)
rod1_A = Line(pivot_point, pivot_point, stroke_color=COLOR_ROD, stroke_width=3.5)
rod2_A = Line(pivot_point, pivot_point, stroke_color=COLOR_ROD, stroke_width=3.5)
bob1_A = Dot(radius=0.14, color=COLOR_CYAN)
bob2_A = Dot(radius=0.18, color=COLOR_NEON_A)
# Setup Visual Elements for System B (Electric Pink)
rod1_B = Line(pivot_point, pivot_point, stroke_color=COLOR_ROD, stroke_width=2.5, stroke_opacity=0.6)
rod2_B = Line(pivot_point, pivot_point, stroke_color=COLOR_ROD, stroke_width=2.5, stroke_opacity=0.6)
bob1_B = Dot(radius=0.12, color="#B5179E", fill_opacity=0.8)
bob2_B = Dot(radius=0.16, color=COLOR_NEON_B)
# Traced paths with glowing effects
trail_A = TracedPath(bob2_A.get_center, stroke_color=COLOR_NEON_A, stroke_width=3.0, dissipating_time=3.5)
trail_B = TracedPath(bob2_B.get_center, stroke_color=COLOR_NEON_B, stroke_width=3.0, dissipating_time=3.5)
# Info Box (Delta theta metric)
delta_label = Text("Perturbation Δθ₀ = 0.001 rad (0.057°)", font="Sans", font_size=18, color=COLOR_ACCENT)
delta_label.to_corner(UL, buff=0.6)
legend_A = Text("● System A (Reference)", font="Sans", font_size=16, color=COLOR_NEON_A)
legend_B = Text("● System B (Perturbed)", font="Sans", font_size=16, color=COLOR_NEON_B)
legend_group = VGroup(legend_A, legend_B).arrange(DOWN, aligned_edge=LEFT, buff=0.15).next_to(delta_label, DOWN, buff=0.25)
self.play(FadeIn(delta_label), FadeIn(legend_group), run_time=1.0)
self.add(trail_A, trail_B, rod1_A, rod2_A, bob1_A, bob2_A, rod1_B, rod2_B, bob1_B, bob2_B)
# Simulation loop state
current_state_A = state_A.copy()
current_state_B = state_B.copy()
dt = 0.02 # Time step
def update_pendulum(mob, *args, **kwargs):
nonlocal current_state_A, current_state_B
# Run sub-steps for higher numerical precision
sub_steps = 4
sub_dt = dt / sub_steps
for _ in range(sub_steps):
current_state_A = rk4_step(current_state_A, sub_dt)
current_state_B = rk4_step(current_state_B, sub_dt)
# Positions A
th1_A, _, th2_A, _ = current_state_A
p1_A = pivot_point + SCALE * np.array([L1 * np.sin(th1_A), -L1 * np.cos(th1_A), 0])
p2_A = p1_A + SCALE * np.array([L2 * np.sin(th2_A), -L2 * np.cos(th2_A), 0])
rod1_A.set_points_by_ends(pivot_point, p1_A)
bob1_A.move_to(p1_A)
rod2_A.set_points_by_ends(p1_A, p2_A)
bob2_A.move_to(p2_A)
# Positions B
th1_B, _, th2_B, _ = current_state_B
p1_B = pivot_point + SCALE * np.array([L1 * np.sin(th1_B), -L1 * np.cos(th1_B), 0])
p2_B = p1_B + SCALE * np.array([L2 * np.sin(th2_B), -L2 * np.cos(th2_B), 0])
rod1_B.set_points_by_ends(pivot_point, p1_B)
bob1_B.move_to(p1_B)
rod2_B.set_points_by_ends(p1_B, p2_B)
bob2_B.move_to(p2_B)
# Attach updater and animate
rod1_A.add_updater(update_pendulum)
# Phase 1: Near-identical synchronized motion (t: 0 to 4s)
sync_text = Text("Initially: Identical trajectory", font="Sans", font_size=20, color=COLOR_TEXT).to_edge(DOWN, buff=0.8)
self.play(FadeIn(sync_text), run_time=0.8)
self.wait(3.5)
# Phase 2: Divergence & Chaos onset (t: 4 to 12s)
chaos_text = Text("Chaos Unfolds: Exponential Divergence (e^λt)", font="Sans", font_size=22, color="#FF4B4B").to_edge(DOWN, buff=0.8)
self.play(Transform(sync_text, chaos_text), run_time=1.0)
self.wait(7.0)
# Phase 3: Total Phase Space Exploration
attractor_text = Text("Strange Attractor in 4D Phase Space", font="Sans", font_size=22, color=COLOR_PIVOT).to_edge(DOWN, buff=0.8)
self.play(Transform(sync_text, attractor_text), run_time=1.0)
self.wait(6.0)
# Clean finish
rod1_A.clear_updaters()
self.wait(1.0)
class DoublePendulumShorts(Scene):
"""
Vertical 9:16 layout optimized for YouTube Shorts, Reels, and Mobile.
Aspect Ratio 9:16, centered composition, high visual impact.
"""
def construct(self):
self.camera.background_color = COLOR_BG
# Vertical Hook Title
hook = Text("CAN 0.001° CAUSE CHAOS?", font="Sans", font_size=34, weight=BOLD, color=COLOR_ACCENT)
hook.to_edge(UP, buff=1.0)
badge = Text("The Double Pendulum Effect", font="Sans", font_size=20, color=COLOR_PIVOT)
badge.next_to(hook, DOWN, buff=0.25)
self.play(FadeIn(hook, scale=1.1), FadeIn(badge), run_time=1.0)
# Pivot point positioned for 9:16 vertical framing
SCALE = 0.82
pivot_point = UP * 1.8
pivot_dot = Dot(pivot_point, radius=0.15, color=COLOR_PIVOT)
self.add(pivot_dot)
# Initial conditions: 0.001 rad difference
init_th1 = np.pi / 2.0
init_th2 = np.pi / 2.0
delta_init = 0.001
state_A = np.array([init_th1, 0.0, init_th2, 0.0], dtype=float)
state_B = np.array([init_th1 + delta_init, 0.0, init_th2, 0.0], dtype=float)
rod1_A = Line(pivot_point, pivot_point, stroke_color=COLOR_ROD, stroke_width=4.0)
rod2_A = Line(pivot_point, pivot_point, stroke_color=COLOR_ROD, stroke_width=4.0)
bob1_A = Dot(radius=0.16, color=COLOR_CYAN)
bob2_A = Dot(radius=0.22, color=COLOR_NEON_A)
rod1_B = Line(pivot_point, pivot_point, stroke_color=COLOR_ROD, stroke_width=3.0, stroke_opacity=0.7)
rod2_B = Line(pivot_point, pivot_point, stroke_color=COLOR_ROD, stroke_width=3.0, stroke_opacity=0.7)
bob1_B = Dot(radius=0.14, color="#B5179E", fill_opacity=0.8)
bob2_B = Dot(radius=0.20, color=COLOR_NEON_B)
trail_A = TracedPath(bob2_A.get_center, stroke_color=COLOR_NEON_A, stroke_width=4.0, dissipating_time=4.0)
trail_B = TracedPath(bob2_B.get_center, stroke_color=COLOR_NEON_B, stroke_width=4.0, dissipating_time=4.0)
self.add(trail_A, trail_B, rod1_A, rod2_A, bob1_A, bob2_A, rod1_B, rod2_B, bob1_B, bob2_B)
current_state_A = state_A.copy()
current_state_B = state_B.copy()
dt = 0.02
def update_pendulum(mob, *args, **kwargs):
nonlocal current_state_A, current_state_B
sub_steps = 4
sub_dt = dt / sub_steps
for _ in range(sub_steps):
current_state_A = rk4_step(current_state_A, sub_dt)
current_state_B = rk4_step(current_state_B, sub_dt)
# System A
th1_A, _, th2_A, _ = current_state_A
p1_A = pivot_point + SCALE * np.array([L1 * np.sin(th1_A), -L1 * np.cos(th1_A), 0])
p2_A = p1_A + SCALE * np.array([L2 * np.sin(th2_A), -L2 * np.cos(th2_A), 0])
rod1_A.set_points_by_ends(pivot_point, p1_A)
bob1_A.move_to(p1_A)
rod2_A.set_points_by_ends(p1_A, p2_A)
bob2_A.move_to(p2_A)
# System B
th1_B, _, th2_B, _ = current_state_B
p1_B = pivot_point + SCALE * np.array([L1 * np.sin(th1_B), -L1 * np.cos(th1_B), 0])
p2_B = p1_B + SCALE * np.array([L2 * np.sin(th2_B), -L2 * np.cos(th2_B), 0])
rod1_B.set_points_by_ends(pivot_point, p1_B)
bob1_B.move_to(p1_B)
rod2_B.set_points_by_ends(p1_B, p2_B)
bob2_B.move_to(p2_B)
rod1_A.add_updater(update_pendulum)
# Bottom Call-to-action
cta = Text("Source code & math proof at:\nmath.praman.blog", font="Sans", font_size=20, color=COLOR_TEXT, line_spacing=1.3)
cta.to_edge(DOWN, buff=1.2)
self.wait(4.0)
self.play(FadeIn(cta, shift=UP*0.3), run_time=1.0)
self.wait(10.0)
rod1_A.clear_updaters()
self.wait(0.5)
Resources & Source Code
Feel free to clone, modify, and render these animations for your own projects, videos, or classrooms!
Leave a Reply