Devlog #81 – Wave 61: Double Pendulum Ensemble, Viscoelastic Fluid & DFT/STFT Visualiser

New Simulations

⚖️ Double Pendulum Ensemble — Quantifying the Butterfly Effect

Why an ensemble?

A single double-pendulum simulation shows chaos visually — the unpredictable motion — but tells you nothing about how fast nearby trajectories diverge. An ensemble of N=30 pendulums, all starting from the same angle except for a tiny offset Δθ in θ₁, makes that divergence rate measurable. It is the numerical equivalent of Lorenz's butterfly: change the initial state by 10⁻⁶ rad and the two trajectories will be completely uncorrelated after ~15 seconds.

Lagrangian equations (m₁=m₂=l₁=l₂=1, g=9.81)

denom = 3 - cos(2θ₁ - 2θ₂)

θ₁'' = [ -3g·sinθ₁ - g·sin(θ₁-2θ₂) - 2sin(θ₁-θ₂)·(θ₂'² + θ₁'²cos(θ₁-θ₂)) ] / denom

θ₂'' = [ 2sin(θ₁-θ₂)·(2θ₁'² + 2g·cosθ₁ + θ₂'²cos(θ₁-θ₂)) ] / denom

Each member is integrated independently with RK4 at Δt = 0.005 s. The reference member (k=0) uses exactly θ₁ = θ₁₀. Member k uses θ₁₀ + k·Δθ. The ensemble spread Δθ₁(t) = max|θ₁(k,t) − θ₁(0,t)| is tracked every frame.

Lyapunov exponent estimation

For a chaotic system, the average rate of growth of infinitesimal separations is δ(t) ~ δ(0) · e^λt, so log δ grows linearly with slope λ (bits/s when using log₂, nats/s for ln). The simulation accumulates log₁₀|Δθ|_max over time, then fits a linear regression through the middle 75% of the history (discarding the initial transient and the saturation plateau). The slope × log₂(10) gives λ in bits/s. Positive λ confirms chaos; for regular initial conditions it converges toward zero.

Colour coding and stats

• Ensemble fan — Each member is drawn from hsl(270°…60°), progressively shifting from violet (reference) to yellow (most offset). The fan width is the immediate visual proxy for separation.
• Coherent count — Members whose |Δθ₁| is below 0.1 rad are counted; when this drops to 1 the system is in fully chaotic regime.
• Regime label — Coherent (all overlap), Diverging (fan spreading), Chaotic (saturated separation).
• Presets — Small spread (Δθ = 10⁻⁴), Large spread (Δθ = 0.2), Crossing angles (θ₂ = 2.5), Ultra-fine (Δθ = 10⁻⁷).

💧 Viscoelastic Fluid — Rheology of Polymer Materials

Viscoelasticity: the middle ground

Most real materials are neither purely elastic (springs that remember their shape) nor purely viscous (fluids that forget it). Polymers, hydrogels, biological tissues, silly putty and paint all show viscoelastic behaviour: they respond elastically to fast deformations (high Deborah number) and viscously to slow ones (low Deborah number). The Deborah number De = τ_R/t_obs is the key dimensionless ratio, where τ_R = η/G₀ is the Maxwell relaxation time.

The three models and their equations

All three models describe a material element subject to stress σ and strain ε:

• Maxwell model (spring + dashpot in series): stress relaxes exponentially at constant strain. G(t) = G₀ · e^−t/τᴿ. Purely viscous at long times: no creep limit — strain grows unbounded. G′(ω) = G₀(ωτ)²/[1+(ωτ)²], G″(ω) = G₀ωτ/[1+(ωτ)²].
• Kelvin-Voigt model (spring + dashpot in parallel): no stress relaxation, but creep saturates. J(t) = (1 − e^−t/τᴿ) / G₀. Purely elastic at high frequencies: G′(ω) = G₀, G″(ω) = ηω.
• Standard Linear Solid (SLS / Zener model): adds a third element (spring in parallel with Maxwell arm). Combines bounded creep and proper stress relaxation. G(t) = G₀[α + (1−α)e^−t/τᴿ] where α is the equilibrium modulus ratio. Most accurate for polymers and biological tissues.

Particle simulation

The left canvas shows a 10×10 grid of mass points connected by Maxwell springs (spring G₀ in series with dashpot η). Each bond tracks its internal viscous extension L_v:

F = G₀_bond × (L − L₀ − Lv)       // spring force
dLv/dt = F / η_bond                 // dashpot relaxation (τR = η/G₀)
overdamped: ζ × v = ΣF_springs      // no inertia; instantaneous velocity

When you drag a particle, nearby bonds stretch (blue → red by stress magnitude) and then slowly relax back after you release. The overall blob behaviour is Maxwell-like: the faster you deform it, the stiffer it feels.

🎵 DFT & STFT Visualiser — From Sinusoids to Spectrograms

The Discrete Fourier Transform

Given N time-domain samples x[n], the DFT computes N complex coefficients:

X[k] = Σₙ x[n] · w[n] · e^{−j2πkn/N}    k = 0, 1, …, N/2

where w[n] is a window function. |X[k]| is the magnitude of the k-th frequency bin, with frequency f_k = k · f_s/N Hz. The visualiser computes this directly (O(N²)) for N ≤ 1024 at 60 fps — fast enough for real-time display.

Window functions and spectral leakage

Non-integer cycles in the DFT window cause energy to "leak" from sharp spectral peaks into neighbouring bins. The choice of window is a trade-off between main-lobe width (frequency resolution) and side-lobe level (leakage):

• Rectangular — Highest resolution, worst leakage (−13 dB first side-lobe). Best when the signal perfectly fills the window.
• Hanning — w[n] = 0.5(1 − cos(2πn/N)). First side-lobe at −31 dB. A standard choice for audio.
• Hamming — w[n] = 0.54 − 0.46cos(2πn/N). First side-lobe at −41 dB but non-zero endpoints → slight discontinuity.
• Blackman — Three-term cosine sum, first side-lobe at −57 dB. Widest main lobe; best for detecting weak tones near strong ones.

STFT spectrogram

The Short-Time Fourier Transform applies the DFT to successive overlapping frames: hop = N − overlap samples apart. The result is a 2D matrix: columns are time frames, rows are frequency bins. The simulation maintains a ring buffer of 460 columns (matching the canvas width) colourised with a viridis-like palette mapping DB value (−80…0) to dark blue → cyan → yellow.

Try the chirp preset (frequency sweeping 100 → 800 Hz over 1 second) and watch the diagonal yellow ridge sweep upward on the spectrogram. Then switch to the multi-tone preset and see the three horizontal lines corresponding to the fundamental and harmonics. Changing the window type while viewing a pure sine with the Rectangular window dramatically shows spectral leakage: frequency sidebands appear and disappear as you cycle through windows.

Signal types

• Sine — Single frequency; one spike in DFT. Ideal for demonstrating window leakage.
• Square — Fourier series 4/π Σ sin((2k−1)2πft)/(2k−1); odd harmonics visible as decreasing peaks.
• Sawtooth — All harmonics at amplitude 1/k; richer harmonic content than square wave.
• Chirp — Linear frequency sweep from f₀ to 8f₀ per second; diagonal stripe on spectrogram.
• AM signal — Amplitude modulation: (1 + 0.8sin(2πf_mt))cos(2πf_ct); two sidebands visible around carrier.
• White noise — Flat spectrum; all frequency bins have equal (random) amplitude.
• Multi-tone — Fundamental f₀ plus harmonics at 3f₀, 5f₀, 7f₀ with decreasing amplitude.

Up Next

Wave 62 will continue expanding the library. Candidates under consideration include fluid dynamics simulations, statistical mechanics visualisations, and additional interactive physics demonstrations. Stay tuned.