Category: Uncategorized

In 1963, meteorologist Edward Lorenz sought to model atmospheric convection using simplified fluid equations. Instead of periodic weather, his computer stumbled upon an infinite, non-repeating geometric object living in 3D phase space: The Lorenz Strange Attractor. In this deep dive, we derive the Lorenz equations from the Navier-Stokes equations, prove phase space volume contraction, and simulate a 36-particle swarm in 3D using Manim.

Figure 1: 3D orbital simulation of a 36-particle ensemble tracing out the fractal butterfly lobes in ℝ³.

1. From Navier-Stokes to Rayleigh-Bénard Convection

Consider a horizontal fluid layer of depth $H$ heated from below and cooled from above, creating a temperature gradient $\Delta T$. Under the Boussinesq approximation (where density variations are neglected except in the buoyancy force), the governing equations are:

where $\mathbf{u} = (u, 0, w)$ is velocity, $\theta$ is temperature deviation, $\nu$ is kinematic viscosity, and $\kappa$ is thermal diffusivity.

2. Fourier Truncation: Lorenz’s 3-Mode Model

Lorenz expanded the stream function $\psi(x, z, t)$ (where $u = -\partial \psi / \partial z$ and $w = \partial \psi / \partial x$) and the temperature perturbation $\theta(x, z, t)$ into spatial Fourier modes, truncating at the lowest non-linear orders:

$$\psi(x, z, t) = \frac{\sqrt{2}\kappa(1 + a^2)}{a} X(t) \sin\left(\frac{\pi a x}{H}\right) \sin\left(\frac{\pi z}{H}\right)$$
$$\theta(x, z, t) = \frac{R_c}{\pi R} \Delta T \left[ \sqrt{2} Y(t) \cos\left(\frac{\pi a x}{H}\right) \sin\left(\frac{\pi z}{H}\right) – Z(t) \sin\left(\frac{2\pi z}{H}\right) \right]$$

Substituting these into the fluid equations and projecting onto each orthogonal Fourier mode yields the celebrated dimensionless Lorenz Equations:

$$\dot{x} = \sigma (y – x)$$
$$\dot{y} = x(\rho – z) – y$$
$$\dot{z} = xy – \beta z$$

3. Mathematical Proof: Zero Volume in Phase Space

Consider a parcel of initial conditions in 3D state space $\mathbb{R}^3$ occupying a volume $V(t)$. The rate of volume expansion or contraction is determined by the divergence of the flow field $\mathbf{F} = [\dot{x}, \dot{y}, \dot{z}]^T$:

$$\nabla \cdot \mathbf{F} = \frac{\partial \dot{x}}{\partial x} + \frac{\partial \dot{y}}{\partial y} + \frac{\partial \dot{z}}{\partial z}$$
$$\nabla \cdot \mathbf{F} = -\sigma – 1 – \beta = -(\sigma + 1 + \beta)$$

For canonical values ($\sigma = 10, \beta = 8/3, \rho = 28$):

$$\nabla \cdot \mathbf{F} = -10 – 1 – \frac{8}{3} = -\frac{41}{3} \approx -13.67 < 0$$

Applying Liouville’s theorem for dynamical systems:

$$\frac{dV}{dt} = \int_{V(t)} (\nabla \cdot \mathbf{F}) \, dV = -\left(\sigma + 1 + \beta\right) V(t)$$
$$V(t) = V(0) \exp\left[-\left(\sigma + 1 + \beta\right)t\right]$$

4. The 36-Particle Swarm Simulation

In our Manim 3D animation, we initialize 36 particles inside a microscopic sphere of radius $r = 0.04$ centered at $(1.0, 1.0, 20.0)$.

As the 3D orbital camera sweeps across the coordinate plane:

• $t \in [0, 4\text{s}]$: The particles travel as a unified luminous filament. The micro-perturbations remain visually imperceptible.
• $t \in [4, 10\text{s}]$: When the swarm approaches the saddle point at the origin $(0,0,0)$, the nonlinear stretching term $xy$ tears the swarm apart. Some particles swing into the left wing while others spiral into the right wing.
• $t > 10\text{s}$: Complete spatial divergence occurs. The ensemble illuminates the entire two-lobed manifold, visually proving the existence of the strange attractor.

5. Full Manim 3D Community Source Code

Below is the complete, self-contained Python script implementing the 3D camera rotation and particle swarm numerical integration:

"""
3D Lorenz Strange Attractor & Fractal Particle Ensemble
Engine: Manim Community Edition (ThreeDScene)
Visual Standard: 3D Orbiting Camera, 60+ Particle Swarm with Velocity Gradient Trails
"""

import numpy as np
from manim import *

# Aesthetic Palette
COLOR_BG = "#08090C"
SIGMA = 10.0
RHO = 28.0
BETA = 8.0 / 3.0

def lorenz_deriv(pos):
    x, y, z = pos
    dx = SIGMA * (y - x)
    dy = x * (RHO - z) - y
    dz = x * y - BETA * z
    return np.array([dx, dy, dz])

def rk4_step_3d(pos, dt):
    k1 = lorenz_deriv(pos)
    k2 = lorenz_deriv(pos + 0.5 * dt * k1)
    k3 = lorenz_deriv(pos + 0.5 * dt * k2)
    k4 = lorenz_deriv(pos + dt * k3)
    return pos + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)

# Velocity to Color Interpolation
def get_velocity_color(speed):
    # Normalized speed range approximately 5 to 60
    t = np.clip((speed - 8.0) / 45.0, 0.0, 1.0)
    # Color ramp: Deep Cyan -> Electric Violet -> Bright Magenta -> Golden Yellow
    c1 = np.array([0x00, 0xF2, 0xFE]) / 255.0  # Cyan
    c2 = np.array([0x79, 0x28, 0xCA]) / 255.0  # Violet
    c3 = np.array([0xFF, 0x00, 0x7F]) / 255.0  # Pink
    c4 = np.array([0xFF, 0xD1, 0x66]) / 255.0  # Gold

    if t < 0.33:
        alpha = t / 0.33
        rgb = (1 - alpha) * c1 + alpha * c2
    elif t < 0.66:
        alpha = (t - 0.33) / 0.33
        rgb = (1 - alpha) * c2 + alpha * c3
    else:
        alpha = (t - 0.66) / 0.34
        rgb = (1 - alpha) * c3 + alpha * c4
    return rgb_to_color(rgb)


class LorenzAttractor3D(ThreeDScene):
    """
    3D Master Scene with dynamic orbital camera and luminous particle swarm.
    """
    def construct(self):
        self.camera.background_color = COLOR_BG

        # 3D Axes (Subtle, minimalist)
        axes = ThreeDAxes(
            x_range=[-30, 30, 10],
            y_range=[-30, 30, 10],
            z_range=[0, 50, 10],
            x_length=7,
            y_length=7,
            z_length=6,
            axis_config={"stroke_color": "#22252e", "stroke_width": 1.5, "stroke_opacity": 0.5}
        )
        # Shift coordinate center to balance the attractor visually
        axes.shift(IN * 1.5 + DOWN * 0.5)

        # Scale and translation helper: transforms Lorenz (x,y,z) to scene coordinates
        def to_scene_coords(pt):
            # Center of the attractor is roughly (0, 0, 25)
            scale = 0.14
            return np.array([pt[0] * scale, pt[1] * scale, (pt[2] - 25.0) * scale])

        # Cinematic Camera Initialization
        self.set_camera_orientation(phi=65 * DEGREES, theta=-55 * DEGREES, zoom=1.15)

        # Fixed 2D Overlay Title & Equations (Rendered in screen space)
        self.begin_ambient_camera_rotation(rate=0.14)

        title = Text("The Lorenz Strange Attractor", font="Sans", font_size=28, color="#F3F4F6")
        subtitle = Text("Fractal Topology & Deterministic Chaos in ℝ³", font="Sans", font_size=16, color="#00F2FE")
        title_group = VGroup(title, subtitle).arrange(DOWN, aligned_edge=LEFT, buff=0.15).to_corner(UL, buff=0.5)

        eq1 = Text("dx/dt = σ(y - x)", font="Monospace", font_size=15, color="#9CA3AF")
        eq2 = Text("dy/dt = x(ρ - z) - y", font="Monospace", font_size=15, color="#9CA3AF")
        eq3 = Text("dz/dt = xy - βz", font="Monospace", font_size=15, color="#9CA3AF")
        eq_group = VGroup(eq1, eq2, eq3).arrange(DOWN, aligned_edge=LEFT, buff=0.12).next_to(title_group, DOWN, buff=0.3)

        self.add_fixed_in_frame_mobjects(title_group, eq_group)
        self.play(FadeIn(title_group, shift=DOWN*0.2), FadeIn(eq_group, shift=DOWN*0.2), run_time=1.2)

        # Ensemble of 35 particles initialized in a tight sphere of radius 0.05 around (1.0, 1.0, 20.0)
        NUM_PARTICLES = 36
        np.random.seed(42)
        base_origin = np.array([1.0, 1.0, 20.0])
        perturbations = np.random.normal(0, 0.04, size=(NUM_PARTICLES, 3))
        
        # Color spectrum for initial swarm
        colors = [
            interpolate_color(ManimColor("#00F2FE"), ManimColor("#FF007F"), i / NUM_PARTICLES)
            for i in range(NUM_PARTICLES)
        ]

        # Create Particle Dots and Traced Trails
        particles = []
        trails = []
        states = []

        for i in range(NUM_PARTICLES):
            init_pos = base_origin + perturbations[i]
            states.append(init_pos)
            scene_pos = to_scene_coords(init_pos)
            
            # Particle dot
            dot = Dot3D(point=scene_pos, radius=0.06, color=colors[i])
            # High-fidelity glowing trail
            trail = TracedPath(
                dot.get_center,
                stroke_color=colors[i],
                stroke_width=2.2,
                dissipating_time=4.0,
                stroke_opacity=0.85
            )
            particles.append(dot)
            trails.append(trail)
            self.add(trail, dot)

        # Numerical Integration Updater
        dt = 0.015
        sub_steps = 4
        sub_dt = dt / sub_steps

        def update_swarm(mob, *args, **kwargs):
            nonlocal states
            for i in range(NUM_PARTICLES):
                cur = states[i]
                for _ in range(sub_steps):
                    cur = rk4_step_3d(cur, sub_dt)
                states[i] = cur
                particles[i].move_to(to_scene_coords(cur))

        # Attach updater to root object
        axes.add_updater(update_swarm)
        self.add(axes)

        # Orbit and observe the fractal butterfly wings emerge
        # Phase 1: Close clustering
        self.wait(4.0)

        # Phase 2: Divergence into both lobes
        phase2_text = Text("Swarm Divergence Across Lobes", font="Sans", font_size=18, color="#FFD166")
        phase2_text.to_edge(DOWN, buff=0.6)
        self.add_fixed_in_frame_mobjects(phase2_text)
        self.play(FadeIn(phase2_text), run_time=0.8)
        self.wait(8.0)

        # Phase 3: Dynamic camera tilt sweep
        self.move_camera(phi=75 * DEGREES, theta=60 * DEGREES, zoom=1.25, run_time=6.0)
        self.wait(6.0)

        # Clean wrap up
        axes.clear_updaters()
        self.stop_ambient_camera_rotation()
        self.wait(1.0)


class LorenzAttractorShorts(ThreeDScene):
    """
    Vertical 9:16 Ultra-Cinematic Scene for YouTube Shorts & Reels.
    Tightly framed, fast visual hook, vertical text layout.
    """
    def construct(self):
        self.camera.background_color = COLOR_BG

        def to_scene_coords(pt):
            scale = 0.13
            return np.array([pt[0] * scale, pt[1] * scale, (pt[2] - 25.0) * scale])

        self.set_camera_orientation(phi=70 * DEGREES, theta=-45 * DEGREES, zoom=1.0)
        self.begin_ambient_camera_rotation(rate=0.18)

        # Top Vertical Hook
        hook = Text("THE BUTTERFLY OF CHAOS", font="Sans", font_size=28, weight=BOLD, color="#FFD166")
        hook.to_edge(UP, buff=1.0)
        sub = Text("35 Parallel Realities in ℝ³", font="Sans", font_size=16, color="#00F2FE")
        sub.next_to(hook, DOWN, buff=0.18)

        cta = Text("Math proof & Python code at:\nmath.praman.blog", font="Sans", font_size=18, color="#FFFFFF", line_spacing=1.3)
        cta.to_edge(DOWN, buff=1.0)

        self.add_fixed_in_frame_mobjects(hook, sub, cta)
        self.play(FadeIn(hook, scale=1.1), FadeIn(sub), run_time=1.0)

        NUM_PARTICLES = 36
        np.random.seed(99)
        base_origin = np.array([0.1, 1.0, 10.0])
        perturbations = np.random.normal(0, 0.05, size=(NUM_PARTICLES, 3))

        colors = [
            interpolate_color(ManimColor("#00F2FE"), ManimColor("#FF007F"), i / NUM_PARTICLES)
            for i in range(NUM_PARTICLES)
        ]

        particles = []
        trails = []
        states = []

        dummy_anchor = Dot(ORIGIN, radius=0).set_opacity(0)
        self.add(dummy_anchor)

        for i in range(NUM_PARTICLES):
            init_pos = base_origin + perturbations[i]
            states.append(init_pos)
            scene_pos = to_scene_coords(init_pos)
            
            dot = Dot3D(point=scene_pos, radius=0.07, color=colors[i])
            trail = TracedPath(
                dot.get_center,
                stroke_color=colors[i],
                stroke_width=2.6,
                dissipating_time=4.5,
                stroke_opacity=0.9
            )
            particles.append(dot)
            trails.append(trail)
            self.add(trail, dot)

        dt = 0.016
        sub_steps = 4
        sub_dt = dt / sub_steps

        def update_swarm(mob, *args, **kwargs):
            nonlocal states
            for i in range(NUM_PARTICLES):
                cur = states[i]
                for _ in range(sub_steps):
                    cur = rk4_step_3d(cur, sub_dt)
                states[i] = cur
                particles[i].move_to(to_scene_coords(cur))

        dummy_anchor.add_updater(update_swarm)
        self.wait(18.0)

        dummy_anchor.clear_updaters()
        self.stop_ambient_camera_rotation()
        self.wait(0.5)

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$$

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.

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