Devlog #88 – Wave 67: Logistic Map, Moiré Pattern, Delaunay Triangulation, Traffic Jam, Quadtree & Tides

New Simulations

📉 Logistic Map — A Window into Universal Chaos

The map and its bifurcations

The logistic map is the simplest 1D system that exhibits the full period-doubling route to chaos:

x_{n+1} = r · x_n · (1 − x_n)      x ∈ [0, 1], r ∈ [0, 4]

Fixed point (period-1 orbit):     x* = 1 − 1/r     stable for r < 3
Period-2 bifurcation at r₁ ≈ 3.000
Period-4 bifurcation at r₂ ≈ 3.449
Period-8 bifurcation at r₃ ≈ 3.544
...
Accumulation point (onset of chaos) at r∞ ≈ 3.56995

To render the bifurcation diagram the simulation iterates the map 300 times (transient) then records the next 300 iterates for each r value, plotting a vertical column of attractor points. The horizontal axis spans r = 2.5 to 4.0 at 1200 pixel resolution. Drag-to-zoom re-renders the selected rectangle at full resolution, revealing self-similar structure inside each chaotic window.

The Feigenbaum constant

Mitchell Feigenbaum discovered in 1975 that the ratio of successive bifurcation intervals converges to a universal constant, independent of the specific map:

Feigenbaum constant:
  δ = lim_{n→∞} (r_n − r_{n-1}) / (r_{n+1} − r_n) ≈ 4.66920160910299...

Convergence:
  (r₁ − r₀) / (r₂ − r₁) = (3.449 − 3.000) / (3.544 − 3.449) = 0.449 / 0.095 ≈ 4.73
  (r₂ − r₁) / (r₃ − r₂) ≈ 4.66  (already very close to the limit)

The same δ appears in the quadratic map, the sine map, and any unimodal map
with a quadratic maximum — a universality class of chaotic systems.

The simulation annotates each period-doubling bifurcation point with a vertical dashed line and computes the current ratio between consecutive intervals, displaying the convergence to δ in a small panel at the side.

Cobweb diagram and Lyapunov exponent

The cobweb diagram (spider-web plot) is a geometric tool for visualising orbit dynamics on the parabola y = rx(1−x). Starting at x₀, draw a vertical line to the parabola at (x₀, x₁), then a horizontal line to the diagonal y = x at (x₁, x₁), then vertical to the parabola at (x₁, x₂), and so on. A stable fixed point produces a spiral converging to the intersection; a period-2 orbit produces a rectangle; chaos produces a dense web.

Lyapunov exponent (finite-time estimate over N iterates):
  λ = (1/N) · Σ_{n=0}^{N-1}  ln |df/dx at x_n|
    = (1/N) · Σ_{n=0}^{N-1}  ln |r(1 − 2x_n)|

  λ < 0 : stable periodic orbit (nearby trajectories converge)
  λ = 0 : bifurcation point or marginal stability
  λ > 0 : chaos (nearby trajectories diverge exponentially)

The Lyapunov exponent plot is drawn directly below the bifurcation diagram, sharing the same r-axis. The chaotic bands (λ > 0, shown in red) and the periodic windows (λ < 0, shown in blue) align exactly with the bifurcation structure above — the most striking being the large period-3 window at r ≈ 3.83, predicted by the Li-Yorke "period-3 implies chaos" theorem.

🌐 Moiré Pattern — Interference, Sampling and Aliasing

Spatial frequency and the moiré period

A line grating with spatial frequency f (lines per pixel) can be described as a periodic binary function g(x) = round(sin(2πfx)). When two gratings with frequencies f₁ and f₂ are multiplied (or equivalently, one is viewed through the other), the product contains sum and difference frequencies:

g₁(x) · g₂(x) = cos(2πf₁x) · cos(2πf₂x)
              = (1/2)[cos(2π(f₁+f₂)x) + cos(2π(f₁−f₂)x)]

Moiré frequency:    f_m = |f₁ − f₂|
Moiré period:       d_m = 1 / f_m    (pixels per moiré band)

Example: f₁ = 0.10 lines/px, f₂ = 0.09 lines/px
  f_m = 0.01 lines/px  →  d_m = 100 px (a wide, visible moiré band)

The simulation renders each grating on a separate off-screen canvas and composites them using globalCompositeOperation = 'multiply'. The live readout shows f₁, f₂, f_m, and d_m updating in real time as you drag the frequency sliders.

Rotated gratings

When two gratings with identical frequency f are rotated by angle θ relative to each other, the moiré pattern has a period determined by the rotation:

Moiré period for equal-frequency gratings at angle θ:
  d_m = d / (2 · sin(θ/2))

  At θ = 5°:  d_m = d / (2 · sin(2.5°)) ≈ 11.5 d   (broad bands)
  At θ = 30°: d_m = d / (2 · sin(15°))  ≈ 1.93 d   (narrow bands)
  At θ = 0°:  d_m → ∞  (parallel gratings: no moiré, uniform field)

A rotation dial lets you sweep from 0° to 90°. At small angles the moiré bands are wide and slowly oscillating; at 45° you get a diamond grid pattern; at 90° the two gratings are orthogonal and produce a square mesh. The moiré period label updates continuously.

Connection to Nyquist aliasing

Moiré in digital images is a direct manifestation of aliasing: when an image sensor samples a fine pattern at a rate below the Nyquist limit, the high spatial frequencies "fold back" to lower frequencies — appearing as coarser patterns that were not in the original scene. The simulation's aliasing demo renders a fine sinusoidal target at frequency f, then downsamples it at rate f_s:

Nyquist criterion: f_s ≥ 2 · f   (no aliasing)

If f_s < 2f (undersampling):
  Alias frequency = |f − round(f / f_s) · f_s|

Example: f = 0.45 cycles/px, f_s = 1 sample/px
  Alias = |0.45 − 0| = 0.45  (no alias, barely inside Nyquist)

  f = 0.55 cycles/px, f_s = 1 sample/px
  Alias = |0.55 − 1| = 0.45  (aliased to 0.45 — same as above!)

The side-by-side display shows the original pattern (fine grid), the sampled result (moiré), and the alias frequency label. This provides the most intuitive explanation of why cameras need anti-aliasing optical low-pass filters in front of their sensors.

△ Delaunay Triangulation — Optimal Triangles and the Voronoi Dual

The Delaunay criterion

A triangulation of a point set is a Delaunay triangulation if and only if for every triangle, the circumscribed circle (circumcircle) contains no other point from the set in its interior. This "empty circumcircle" criterion is equivalent to maximising the minimum angle across all triangles in the triangulation — making Delaunay triangulations the most "equilateral" possible for a given point set.

Circumcircle of triangle (A, B, C):
  Centre D can be found by intersecting perpendicular bisectors of AB and BC.

  In practice (robust computation):
    ax, ay = A.x − C.x,  A.y − C.y
    bx, by = B.x − C.x,  B.y − C.y
    D_denom = 2 · (ax·by − ay·bx)
    D.x = (by·(ax²+ay²) − ay·(bx²+by²)) / D_denom + C.x
    D.y = (ax·(bx²+by²) − bx·(ax²+ay²)) / D_denom + C.y
    R   = distance(D, A)

The hover inspector draws the circumcircle of the triangle under the cursor, shading it red if it (erroneously) contains another point — a useful debugging mode that verifies the algorithm's correctness in real time.

Bowyer-Watson incremental algorithm

The simulation uses the Bowyer-Watson algorithm, which incrementally inserts points one at a time:

Initialise with a "super-triangle" that contains all points.

For each new point P:
  1. Find all "bad" triangles: triangles whose circumcircle contains P.
  2. Find the boundary polygon of the "bad" triangles:
       edges not shared between two bad triangles.
  3. Delete all bad triangles.
  4. Create new triangles from P to each boundary edge.
  5. Re-check and fix any remaining circumcircle violations (edge flips).

After all points inserted:
  Remove all triangles sharing a vertex with the super-triangle.

Complexity: O(n log n) expected for uniformly random point sets;
             O(n²) worst case (degenerate configurations).

Each insertion step is animated: bad triangles flash red, the cavity polygon is drawn in orange, then the new triangles fan in from the inserted point. A step-by-step mode lets you advance one point at a time to trace the algorithm.

Voronoi diagram as the dual

The Voronoi diagram is the geometric dual of the Delaunay triangulation. For each Delaunay triangle, place a vertex at its circumcenter; connect circumcenters of adjacent triangles (sharing an edge). The result is the Voronoi diagram: each cell contains all points closer to one input point than to any other.

Duality:
  Delaunay edge (A,B) ↔ Voronoi edge separating cells A and B
  Delaunay vertex P   ↔ Voronoi cell of P
  Delaunay triangle   ↔ Voronoi vertex (circumcenter)

The Voronoi cell of point P is the convex polygon:
  { x : |x − P| ≤ |x − Q| for all other input points Q }

Toggle the "Show Voronoi" button to overlay the dual graph. Both structures are drawn simultaneously in different colours, making the duality relationship visually obvious: every Delaunay edge is crossed perpendicularly by the corresponding Voronoi edge.

🚗 Traffic Jam — How Phantom Jams Emerge from Rules

The Nagel-Schreckenberg model

The Nagel-Schreckenberg (NaSch) cellular automaton (1992) models traffic on a one-lane ring road divided into L cells. Each cell is either empty or occupied by one car with integer velocity v ∈ {0, 1, 2, …, v_max}. One time step applies four rules simultaneously to all cars:

Let d = gap ahead (number of empty cells before next car)
    v = current velocity
    v_max = maximum velocity (default 5)
    p = randomisation probability (default 0.3)

Rule 1 — Acceleration:
    if v < v_max:  v ← v + 1

Rule 2 — Braking (collision avoidance):
    if v ≥ d:  v ← d − 1

Rule 3 — Randomisation (human reaction noise):
    if v > 0 and random() < p:  v ← v − 1

Rule 4 — Movement:
    car advances v cells forward

The randomisation step is essential: without it (p = 0) traffic flows freely at any density. With p > 0, random braking events propagate backwards as stop-and-go waves — the "phantom jam" that drivers experience on motorways with no visible cause.

Space-time diagram and the fundamental diagram

The space-time diagram plots car positions (horizontal axis = road position, vertical axis = time, scrolling downward). Each car is a dot; jams appear as dark diagonal stripes slanting to the left — indicating the jam travels backwards at a roughly constant speed. The jam speed in the NaSch model is typically −1 cell/step (one cell backward per time step) regardless of traffic density.

Traffic flow variables:
  ρ = density = N / L  (cars per cell, 0 ≤ ρ ≤ 1)
  q = flow    = ρ · ⟨v⟩  (cars passing a point per time step)

Fundamental diagram (q vs ρ) for NaSch with v_max = 5, p = 0.3:
  Free-flow branch  (ρ < ρ_crit ≈ 0.06):  q ≈ ρ · v_max
  Congested branch  (ρ > ρ_crit):          q decreases with ρ
  Maximum flow at  ρ_crit ≈ 0.06,  q_max ≈ 0.32 cars/step/cell

The live fundamental diagram is plotted as a scatter point (current density, current flow) on the q-ρ plane, which traces out the characteristic inverted-V shape as density is varied with the slider. The hysteresis between the free-flow and congested branches — familiar from empirical motorway data — is reproduced qualitatively by the model.

Phase transition and the faster-is-slower effect

At critical density the model undergoes a first-order phase transition from free flow to congested flow. Above ρ_crit, increasing v_max can paradoxically decrease throughput: faster individual speeds mean larger safe-following distances are required, reducing density, which reduces flow. This "faster-is-slower" effect is also seen in pedestrian crowd models (Helbing) and connects to the broader physics of driven dissipative systems near jamming transitions.

🌲 Quadtree — Spatial Partitioning and Collision Broadphase

Quadtree structure and insertion

A quadtree recursively subdivides a 2D region into four equal quadrants. Each internal node represents a rectangle; each leaf either holds one point (or is empty). Insertion of a new point P:

function insert(node, P):
    if node is leaf and node is empty:
        store P in node
        return
    if node is leaf and node holds point Q:
        subdivide node into 4 children
        reinsert Q into appropriate child
        insert P into appropriate child
        return
    // node is internal:
    insert P into the child quadrant that contains P

Depth for n uniformly random points: O(log n) expected
Worst case: n collinear points → O(n) depth (degenerate)

The canvas visualisation draws every rectangle boundary in real time. As points are added one by one, new subdivisions appear with a brief flash animation. The tree depth and leaf count are shown in the stats panel.

Range query with branch pruning

A range query asks: "which points lie inside this rectangle R?" Naïve search checks all n points in O(n). The quadtree prunes branches whose bounding rectangle does not overlap R:

function range_query(node, R):
    if node.bounds does not overlap R:
        return []    // prune: skip entire subtree
    if node is leaf:
        return [P for P in node if P in R]
    return range_query(NW) + range_query(NE) +
           range_query(SW) + range_query(SE)

Complexity: O(√n + k) for uniformly random points,
            where k = number of results returned.

The interactive range query lets you draw a query rectangle by clicking and dragging. Pruned subtrees are immediately shaded grey; active subtrees glow. A counter shows "nodes visited" vs "total nodes" to quantify the pruning benefit.

Collision detection broadphase

The collision-detection panel simulates n moving circles. The naïve approach tests all n(n−1)/2 pairs each frame. The quadtree broadphase rebuilds the tree each frame (O(n log n)) and only tests circles in the same or adjacent cells as broadphase candidates:

Naïve:     O(n²) checks per frame
Quadtree:  O(n log n) build + O(n·k) checks
           where k = average neighbours per cell ≈ constant

Speedup factor at n = 500:
  Naïve: 124,750 pair checks per frame
  Quadtree: ~2,500 broadphase candidates (≈ 50× fewer checks)

The live counter shows both check counts side by side, updating each frame. The speedup ratio is plotted against n as a secondary graph, confirming the asymptotic O(n²) vs O(n log n) divergence.

🌊 Tides — The Gradient of Gravity

Tidal force: the gradient of gravity

Tidal forces are not gravity itself, but the difference in gravitational acceleration across a finite-size body. For the Earth-Moon system, the tidal acceleration at a point displaced Δr from Earth's centre toward the Moon (distance D) is:

Gravitational acceleration from Moon at Earth's centre:
  a₀ = G·M_Moon / D²

At a surface point displaced Δr toward Moon:
  a = G·M_Moon / (D − Δr)²
  ≈ G·M_Moon / D² · (1 + 2Δr/D)   (Taylor expansion, Δr ≪ D)

Tidal acceleration (differential):
  Δa = a − a₀ = 2·G·M_Moon·Δr / D³   ∝ 1/D³

This 1/D³ dependence means tidal forces fall off faster than gravity (1/D²),
so the Moon's tidal effect on Earth exceeds the Sun's despite the Sun's
greater gravitational pull:
  (Tidal force)_Moon / (Tidal force)_Sun ≈ 2.2

The vector field display draws tidal acceleration arrows at grid points around Earth's surface. The arrows form the characteristic "two-lobe" pattern: outward stretching toward and away from the Moon, inward compression perpendicular to the Earth-Moon line. This gives two tidal bulges and explains the approximately twice-daily high-tide cycle.

Spring and neap tides

The Moon and Sun each produce their own tidal force. Their relative alignment determines the type of tide:

Syzygy (Moon–Earth–Sun or Sun–Earth–Moon aligned):
  Tidal forces add constructively → spring tides
  Occurs at new moon (conjunction) and full moon (opposition)
  Spring tide amplitude ≈ 1 + 0.46 = 1.46 × average

Quadrature (Moon 90° from Sun as seen from Earth):
  Tidal forces partially cancel → neap tides
  Occurs at first quarter and last quarter
  Neap tide amplitude ≈ 1 − 0.46 = 0.54 × average

  (where 0.46 is the Sun/Moon tidal ratio from mass and distance)

The Moon-phase dial sweeps through the 29.5-day synodic month. As it rotates, the vector field superposition updates in real time, and the marigraph below records the resulting tidal range. The spring-tide peaks align with the new and full moon positions; neap troughs fall at the quarters.

Tidal locking

The Moon always presents the same face to Earth — a consequence of tidal locking. Earth's gravity raises a tidal bulge in the Moon's rock. Early in the Moon's history, when it rotated faster, this bulge was dragged ahead of the Earth-Moon line by rotational inertia. The gravitational torque from Earth on this misaligned bulge decelerated the Moon's rotation until its rotation period equalled its orbital period. The simulation shows an animation of this process with the tidal bulge angle and the decelerating torque labelled.

Tidal locking timescale:
  t_lock ∝ a⁶ · ω₀ · Q / (M_planet · R_moon⁵)

  where a = orbital semi-major axis, Q = tidal quality factor,
        ω₀ = initial spin rate, R_moon = radius of tidally locked body

Moon-Earth: t_lock ≈ a few hundred million years (already locked)
Earth-Sun:  t_lock ≈ 50 billion years (will not lock within Sun's lifetime)

Up Next

Wave 68 continues the momentum with six more simulations: Barnes-Hut O(n log n) N-body gravity, marching squares iso-contouring with metaballs, Bridson Poisson disk blue-noise sampling, the CHSH Bell inequality experiment reaching the Tsirelson bound, the Chua circuit double-scroll attractor, and a Helbing social-force pedestrian crowd model with lane formation. That wave and its devlog are already live — see Devlog #89.