Task 6 of 6

Payoff: RK4 Drifts, Verlet Does Not

Eighty-one orbits at dt = 0.5, every trajectory released from x = 0 with v = 1. That start makes the bookkeeping free: E₀ = ½, so the energy as a multiple of its starting value is exactly x² + v².

The launch is the interesting part. Thread x integrates this.thread.x + 1 steps — one thread takes a single step, the last takes all 1,024, and between them the 1,024 threads trace the whole energy history of the run. Every prefix computed independently, from scratch, at the same time. Nothing is shared, so nothing has to be sequenced.

RK4 is the most accurate thing you have written: after eighty-one orbits its phase is 0.24 radians behind the truth, where Verlet's has slipped nearly a whole orbit. Watch what its energy does anyway — and remember that midpoint, second order and perfectly respectable, ends this same run with about eight million times the energy it started with.

the accurate method leaks; the symplectic one ripples and stays

Goal: finish the velocity-Verlet energy kernel, then log each method's final energy and the smallest value Verlet ever reaches.

Requirements

Thread x runs this.thread.x + 1 steps — the loop bound is the thread index
Start every trajectory at x = 0, v = 1 and step with velocity Verlet
Return x * x + v * v — the energy as a multiple of where it started
console.log both final energies and the minimum of the Verlet curve

Hint 1 — the trip count

for (let s = 0; s < this.thread.x + 1; s++). gpu.js compiles a bound it cannot fold at build time into a loopMaxIterations loop with an early exit, which is why the kernel declares loopMaxIterations: 1024.

Hint 2 — why x² + v² is the energy ratio

E = (x² + v²)/2 and every trajectory starts at (0, 1), so E₀ = ½ and E/E₀ = x² + v². A value of 1 means the energy is exactly where it began.

Hint 3 — what to look for in JavaScript

Scan for the smallest value on each curve, and read the last one:

let low = Infinity;
for (let i = 0; i < verletCurve.length; i++) {
  if (verletCurve[i] < low) low = verletCurve[i];
}

RK4's minimum is its final value — it never goes back up.

Same idea elsewhere

Two lessons travel. The launch shape — a thread per prefix, wildly unequal work, the wavefront running until its slowest lane finishes — is the load-imbalance problem every CUDA, ROCm, WebGPU and Metal programmer eventually has to lay out differently. And the result is why orbital-mechanics and molecular-dynamics codes on those platforms pick symplectic integrators: over a long run you are not choosing the method with the smallest error per step, you are choosing the one whose error does not have a direction.

All tasks in ODE Integrators

This page is an interactive exercise — the editor, the GPU runner and your saved progress need JavaScript. The text above is the full brief.