# Payoff: RK4 Drifts, Verlet Does Not

*Task 6 of 6 · [ODE Integrators](https://gpu.rocks/learn/ode-integrators-62f4a3ff.md) · GPU.js Learn*

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.

## Figures

- **the accurate method leaks; the symplectic one ripples and stays**

## Goal

**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:

```js
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.

## Starter code

```js
// 1,024 steps of dt = 0.5 — about eighty-one orbits.
// Thread x integrates x + 1 of them, so the 1,024 threads together
// trace the whole energy history, every prefix computed from scratch.
const gpu = new GPU({ mode });
const DT = 0.5;

// RK4, exactly as you wrote it two tasks ago, reporting energy instead of error.
const rk4Energy = gpu.createKernel(function (dt) {
  let x = 0;
  let v = 1;
  for (let s = 0; s < this.thread.x + 1; s++) {
    const half = dt / 2;
    const k1x = v;
    const k1v = -x;
    const k2x = v + k1v * half;
    const k2v = -(x + k1x * half);
    const k3x = v + k2v * half;
    const k3v = -(x + k2x * half);
    const k4x = v + k3v * dt;
    const k4v = -(x + k3x * dt);
    x = x + (dt / 6) * (k1x + 2 * k2x + 2 * k3x + k4x);
    v = v + (dt / 6) * (k1v + 2 * k2v + 2 * k3v + k4v);
  }
  return x * x + v * v;
}, { output: [1024], loopMaxIterations: 1024 });

const verletEnergy = gpu.createKernel(function (dt) {
  let x = 0;
  let v = 1;
  // TODO: run this.thread.x + 1 velocity-Verlet steps and return
  // the energy as a multiple of its starting value: x * x + v * v.
  return 1;
}, { output: [1024], loopMaxIterations: 1024 });

const rk4Curve = await rk4Energy(DT);
const verletCurve = await verletEnergy(DT);

// TODO: log both final energies, and the smallest value the Verlet curve reaches.
```

---

Interactive version: https://gpu.rocks/learn/ode-integrators-62f4a3ff/6

[Previous task](https://gpu.rocks/learn/ode-integrators-62f4a3ff/5.md)