Devlog #93 – Wave 72: Kalman Filter, Backpropagation, 3D Pendulum, Eutrophication, Acoustic Lens & Planetary Rings

New Simulations

Kalman Filter: Optimal State Estimation

The Kalman filter (1960) is one of the most successful algorithms in applied mathematics. It optimally combines a dynamic model's prediction with noisy sensor measurements, producing the minimum-variance unbiased state estimate — and crucially, it tracks exactly how uncertain it is via the covariance matrix P.

The 4D State and Full Matrix Operations

The state vector is [x, y, v_x, v_y] — position and velocity in 2D. The 4×4 covariance matrix P tracks uncertainties and correlations between all state components. Each predict step computes:

x̂⁻ = F·x̂          (state prediction via transition matrix F)
P⁻  = F·P·Fᵀ + Q   (covariance prediction, Q = process noise)

Each update step (on receiving a 2D position observation via H = [I₂|0₂]):

K  = P⁻·Hᵀ·(H·P⁻·Hᵀ + R)⁻¹  (Kalman gain)
x̂  = x̂⁻ + K·(z − H·x̂⁻)       (state update)
P  = (I − K·H)·P⁻·(I − K·H)ᵀ + K·R·Kᵀ  (Joseph form)

Joseph Form for Numerical Stability

The standard update P = (I−KH)P⁻ suffers from floating-point asymmetry — P can lose positive-definiteness after many iterations. The Joseph form P = (I−KH)P⁻(I−KH)ᵀ + KRKᵀ is algebraically equivalent but numerically stable: it guarantees P remains symmetric positive definite regardless of floating-point rounding.

Covariance Ellipse Visualisation

The 2D positional uncertainty is visualised as an ellipse computed from the eigendecomposition of the 2×2 positional submatrix of P. The semi-axes are √(λ₁) and √(λ₂) (standard deviations) and the orientation follows the eigenvectors. The ellipse visually swells during the predict step and shrinks when an observation arrives — a beautiful demonstration of the filter's information cycle.

• 4 trajectories: straight line, circular orbit, random walk, evasive manoeuvre
• RMSE comparison: Kalman estimate vs raw measurements vs model prediction only
• Process noise Q and measurement noise R sliders
• Toggle observation frequency to see covariance grow during gaps

Backpropagation: Gradient Flow Made Visible

Backpropagation is just the chain rule applied efficiently to a computation graph — but understanding it viscerally requires seeing the gradient pulses flow backward through each layer. This simulation makes that flow explicit.

Manual Chain Rule Implementation

The network architecture is configurable: 2 inputs → H hidden units → H hidden units → 2 (classification) or 1 (regression) outputs. Each forward pass computes activations layer by layer. Each backward pass computes δ errors: output layer gets δ_L = ∂L/∂z_L, then propagates: δ_l = (Wᵀ_{l+1} · δ_{l+1}) ⊙ σ'(z_l). Weight gradients are ∂L/∂W_l = δ_l · aᵀ_{l-1}.

Visual Gradient Animation

During backprop, coloured pulses travel right-to-left through the network diagram. Edge width and colour intensity represent |∂L/∂w| for each weight. Large gradients glow bright; near-zero gradients (vanishing gradient) show as dim lines. With ReLU, dead neurons (always outputting 0) are marked explicitly — their outgoing gradients are permanently zero.

• Activations: ReLU, tanh, sigmoid; output: softmax+CE or linear+MSE
• Decision boundary heatmap updates every N training steps
• Loss strip plot: training loss over time
• Step-through single-sample mode for learning the algorithm

3D Pendulum: Lagrange Multiplier Constraint

A spherical pendulum — a mass on a rigid rod, free to swing in all directions — is simple to describe physically but awkward to simulate numerically. The constraint is that the distance from pivot to mass is always exactly L. Rather than switching to spherical coordinates (which introduce gimbal-lock singularities), the simulation works in Cartesian coordinates and enforces the constraint via a Lagrange multiplier.

Cartesian RK4 with Constraint Projection

The equations of motion in Cartesian (x,y,z) coordinates are: m·r̈ = F_gravity + λ·r̂ where λ is the constraint force magnitude (tension). At each RK4 step, after advancing r and v freely, the positions are projected back to the sphere and velocities to the tangent plane: λ = (−|v|² − r·g_vec) / L², then apply the correction. This ensures exact constraint satisfaction without coordinate singularities.

Three Physical Modes

• Planar rosette: initial conditions in a single plane — traces a Lissajous-like rosette pattern when viewed from above due to the nonlinear period
• Conical pendulum: mass orbits at fixed angle θ with period T = 2π√(L cosθ / g) — shorter than the simple pendulum period due to the reduced effective gravity
• Foucault pendulum: adds Coriolis force F_c = −2mΩ × v with Ω = Ω_Earth·sin(latitude)·ẑ, causing the plane of oscillation to precess at rate Ω·sin(lat)

• 3D projection with mouse drag-to-rotate
• Latitude slider for Foucault mode (0° = equator, 90° = pole)
• Trace length and colour gradient by time

Eutrophication: The Lake Tipping Point

Eutrophication — nutrient enrichment leading to algal bloom and oxygen depletion — is one of the most important examples of a catastrophic tipping point in environmental science. Once a lake crosses the threshold, it snaps to a turbid, low-oxygen state that can persist for decades even after the phosphorus input is reduced.

The Five-Species ODE System

The simulation integrates five coupled ODEs representing P (dissolved phosphorus), A (algae biomass), O (dissolved oxygen), F (fish biomass), and D (detritus). The key nonlinearity is the phosphorus recycling term: under anoxic conditions (O below threshold), sediment releases stored phosphorus — a positive feedback that drives the system to an alternate stable state. The RK4 integrator handles the stiff transitions reliably.

Hysteresis and the Phase Plot

The phase plot traces (phosphorus loading, algae biomass) over time. As loading increases slowly, the system stays in the clear-water state until it abruptly jumps to the bloom state at the upper saddle-node bifurcation point. As loading is then reduced, the system remains in the bloom state all the way down to the lower saddle-node — this is hysteresis. The area enclosed by the hysteresis loop in the phase plot represents the "recovery debt" — how much extra effort is needed to restore the lake compared to preventing the bloom in the first place.

• Tank visualisation: water colour gradient from clear blue to murky green, bloom surface mat, fish icons, sediment darkness
• Phosphorus loading slider with slow-increase and slow-decrease sweep modes
• Real-time oxygen level and fish population readouts

Acoustic Lens: FDTD Wave Simulation

Acoustic focusing — concentrating sound waves to a point — underlies medical ultrasound, sonar, and architectural acoustics. This simulation solves the 2D wave equation directly using the Finite Difference Time Domain (FDTD) method, showing the wavefronts bending in real time.

Leapfrog FDTD and Mur Boundaries

The 2D pressure wave equation ∂²p/∂t² = c²(∂²p/∂x² + ∂²p/∂y²) is discretised on a 240×144 grid using the leapfrog scheme: pressure at time n+1 depends on times n and n-1 and the spatial Laplacian at n. The CFL stability condition requires c·dt/dx ≤ 1/√2. The first-order Mur absorbing boundary condition eliminates reflections from the grid edges by matching the outgoing wave speed: p_edge(t+dt) = p_{edge+1}(t) + (c·dt−dx)/(c·dt+dx)·(p_{edge+1}(t+dt) − p_edge(t)).

Four Lens Types

• Convex lens: higher sound speed in the lens region (lower density material) converges the wavefront
• Concave lens: same material but concave geometry diverges the beam
• Fresnel lens: zone-plate design — alternating high/low speed concentric rings that add in phase at the focal point
• GRIN (parabolic): speed varies as c(r) = c₀ / (1 + (r/r₀)²) — a graded index lens with no discrete boundaries

• Divergent heatmap: blue/red pressure field with additive glow
• Snell's law ray overlay: geometric optics approximation overlaid on wave simulation
• EMA pressure at predicted focal point validates focal length against paraxial formula

Planetary Rings: Keplerian Particles and Resonance Gaps

Saturn's rings are one of the solar system's most photogenic structures — and their structure is controlled by orbital mechanics. The Cassini Division, a prominent gap between the A and B rings, is emptied by the 2:1 orbital resonance with the moon Mimas.

Keplerian Particles and Resonance Kicks

Up to 10,000 test particles are placed on Keplerian orbits with semi-major axis a distributed as n ∝ a^(−3/2) (flat surface density). Each particle receives a gravitational kick from a Mimas-like moon at each conjunction. At the 2:1 resonance radius (where the particle completes exactly 2 orbits per moon orbit), these kicks always occur at the same orbital phase — they accumulate coherently and pump up eccentricity until the orbit becomes crossing, after which the particle is quickly scattered out. After thousands of orbits, the resonance radius shows a density gap matching the Cassini Division.

Roche Limit and Shepherd Moons

The Roche limit R_Roche = a_planet · (2ρ_planet/ρ_satellite)^(1/3) marks where tidal forces overcome self-gravity — moons inside it are shredded into ring material. The simulation also supports shepherd moon pairs: two small moons straddling a narrow ringlet exchange angular momentum with ring particles, keeping the ringlet sharply confined — the mechanism behind Saturn's F-ring structure.

• N slider: 100 to 10,000 particles
• Moon mass and orbital radius adjustable — watch the gap form and move
• Colour by eccentricity to visualise resonance excitation
• Toggle shepherd moon pair for narrow ringlet confinement

Up Next

Wave 73 is a WebGL/GLSL special: all five simulations use custom shaders, raymarching, or GPGPU techniques. Coming up: Mandelbulb 3D fractal raymarching with distance estimator, Menger sponge SDF with Hausdorff dimension, Morris-Thorne wormhole Flamm's paraboloid embedding, Gielis supershape spherical product, and 250,000-particle GPGPU simulation with ping-pong render targets.