The Two Shear Transforms
Every timestep decomposes into two sequential maps. Each is a linear shear — a parallelogram-preserving tilt in phase space. Their composition is the full symplectic step.
T₁ — Kick (shear in velocity space)
─────────────────────────────────────
Compute acceleration from position:
aₓ = −GM · x / r³
aᵧ = −GM · y / r³ (r = √(x²+y²))
Update velocity first:
vₓ ← vₓ + aₓ(x, y) · Δt
vᵧ ← vᵧ + aᵧ(x, y) · Δt
T₂ — Drift (shear in position space)
─────────────────────────────────────
Update position using new velocity:
x ← x + vₓ_new · Δt
y ← y + vᵧ_new · Δt
This is the crucial "line flip":
T₂ consumes T₁'s output, not old v.
Jacobian Proof: det = 1
Each map is area-preserving because its Jacobian has unit determinant.
T₁ (force ⊥ to velocity space):
∂(vₓ',vᵧ')/∂(vₓ,vᵧ) = I₂ → det = 1 ✓
(force depends only on x,y not on v)
T₂ (drift ⊥ to position space):
∂(x',y')/∂(x,y) = I₂ → det = 1 ✓
(shift depends only on v not on x,y)
Full map T = T₂ ∘ T₁:
det(J_T) = det(J_T₂)·det(J_T₁) = 1·1 = 1
⟹ 4-D phase-space volume preserved
at every discrete step, for all time.
Naive Euler — Why It Fails
The red particle uses simultaneous updates (old v → x, old a → v). Both branches use the old state independently:
Naive Euler (simultaneous, unstable)
─────────────────────────────────────
x_new = x + v_old · Δt ← old v
v_new = v + a(x_old) · Δt← old a
The full Jacobian det > 1 (area grows).
Phase-space patch expands each step →
energy drifts monotonically outward →
orbit spirals away.
The cyan particle uses the symplectic form — its orbital radius stays bounded indefinitely. Watch the energy chart: symplectic oscillates, naive drifts.
Poincaré Invariants
For forces depending only on position, there is a hierarchy of conserved geometric quantities — the Poincaré integral invariants:
Phase space: (x, vₓ, y, vᵧ) — 4 dimensions
I₁ = ∮ (vₓdx + vᵧdy)
[oriented 2-D area on x–vₓ plane]
I₂ = ∬ dvₓ∧dx + dvᵧ∧dy
[sum of all oriented 2-D sub-areas]
I₃ = ∭∭ dx∧dvₓ∧dy∧dvᵧ
[full 4-D phase-space volume]
All three remain constant under T = T₂∘T₁.
This is Liouville's theorem, discretized.
Energy oscillates in a bounded band rather than drifting — it is NOT conserved exactly, but the geometric invariants that prevent secular drift ARE preserved.
JavaScript Float64 Precision
All particle state is stored in Float64Array — IEEE 754 double precision — 64 bits, 52-bit mantissa:
Machine epsilon: ε ≈ 2.22 × 10⁻¹⁶
Significant digits: ~15–16
In AU/yr units (r ~ 1–6, v ~ 1–10):
Rounding error per step: ~10⁻¹⁵ AU
After N steps (random walk): O(√N · ε)
At N = 10⁶ steps: ~10⁻¹³ AU (sub-atomic)
Symplectic structure keeps errors bounded:
no secular trend → no catastrophic loss
of significant digits over long runs.
Float64Array vs. plain Array:
- Guaranteed 64-bit per element
- Cache-coherent, faster iteration
- No boxing overhead per value
Zoomable Canvas — Affine Transform
World coordinates (AU) map to screen pixels via a camera affine transform updated each frame:
screen_x = W/2 + (world_x − cam_x) · scale
screen_y = H/2 − (world_y − cam_y) · scale
(y-axis flipped: screen grows downward)
Zoom (scroll wheel):
scale ← scale · factor^(−Δwheel)
Pivot at mouse: adjust cam so the
world point under cursor is fixed.
Pan (drag):
cam_x −= Δmouse_x / scale
cam_y += Δmouse_y / scale
Rendering per frame:
1. requestAnimationFrame callback
2. Run stepsPerFrame integration steps
3. Clear canvas, draw grid + stars
4. Batch-render particles (beginPath loop)
5. Draw glowing Sun, Jupiter, refs
6. Update energy deviation chart
Wisdom–Holman Splitting (1991)
The Leapfrog is a stepping stone. The 1991 breakthrough replaced the straight-line Drift (T₂) with an exact Kepler arc — known analytically. Jupiter then delivers a velocity kick between arcs:
T_Kepler: advance along exact ellipse
→ no approximation for Sun's gravity
→ timestep can be much larger
T_kick: apply planetary perturbations
→ small velocity correction
The split is still symplectic → inherits
all Poincaré invariant preservation.
Result: millions of years of stability.
This is how Kirkwood gaps emerge in
simulations that run for ~10⁷ yr.
The gaps appear at orbital resonances with Jupiter (3:1 at 2.50 AU, 5:2 at 2.82, 7:3 at 2.95, 2:1 at 3.27). Labeled on the AU reference circles in the viewport.