# Velocity Verlet: the Symplectic Step

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

**Velocity Verlet** reads like Euler with one extra term. Move the
position with a quadratic instead of a straight line, then update the velocity using the
*average* of the accelerations at the two ends of the step:

```js
x ← x + v·dt + ½·a·dt²
v ← v + ½·(a + a_new)·dt        where a_new = −x_new
```

It is second order — the same exponent midpoint gave you, so the instrument you have
been using *cannot tell them apart* by exponent alone. (It is about four times
more accurate than midpoint at the same step size, and a real implementation carries
`a_new` into the next step as its `a`, so it costs one force
evaluation per step against midpoint's two. Here `a = −x` is free, so the code
below just recomputes it.)

What it has that midpoint does not is invisible in a single step and decisive over a
million. The Verlet step is an exactly area-preserving map of the `(x, v)`
plane — **symplectic** — and for this oscillator it conserves
`(1 − dt²/4)·x² + v²` exactly, forever, in a way no accumulation of steps can
erode. Measure the order here. The next task is where that sentence starts to matter.

## Goal

**Goal:** write the velocity Verlet step. The ratio should be ≈4
again — and the errors themselves about four times smaller than midpoint's.

## Requirements

- Position first, with the quadratic term: `x += v·dt + 0.5·a·dt·dt` where `a = −x` at the start of the step
- Then recompute the acceleration at the NEW position: `aNext = −x`
- Velocity from the average of the two: `v += 0.5·(a + aNext)·dt`

## Hint 1 — order of operations

`a` is captured before `x` moves; `aNext` is
read after. The velocity update needs both, so it has to come last.

## Hint 2 — the extra term

`0.5 * a * dt * dt` is the constant-acceleration formula from
first-year mechanics. Dropping it leaves a scheme that still looks plausible and is
back to first order.

## Hint 3 — the whole step

```js
const a = -x;
x = x + v * dt + 0.5 * a * dt * dt;
const aNext = -x;
v = v + 0.5 * (a + aNext) * dt;
```

## Same idea elsewhere

Velocity Verlet, not RK4, is what every production molecular-dynamics code on a GPU
actually ships — GROMACS, LAMMPS, HOOMD-blue and OpenMM all step with Verlet or leapfrog on
CUDA, HIP and Metal. It is not because they cannot afford RK4's four evaluations. It is
because a run of a hundred million steps is judged on whether the energy stayed put, and
that is a property of the *shape* of the update, not of its order.

## Starter code

```js
// Same instrument again. Same order as midpoint — and something else.
const gpu = new GPU({ mode });

const errorOf = gpu.createKernel(function (startX, startV, levelSteps) {
  const n = levelSteps[this.thread.y];
  const dt = this.constants.tEnd / n;

  let x = startX[this.thread.x];
  let v = startV[this.thread.x];
  for (let s = 0; s < n; s++) {
    // TODO: one velocity-Verlet step.
    //   x moves with v·dt AND the ½·a·dt² term,
    //   then v updates on the AVERAGE of a and the new a.
    const a = -x;
    x = x + v * dt;
    v = v + a * dt;
  }

  const exact = startX[this.thread.x] * Math.cos(this.constants.tEnd)
              + startV[this.thread.x] * Math.sin(this.constants.tEnd);
  return Math.abs(x - exact);
}, {
  output: [256, 4],
  constants: { tEnd: 5 },
  loopMaxIterations: 200,
});

const errors = await errorOf(startX, startV, levelSteps);

// average the trajectories at each step size, then see what each halving bought
const means = [];
for (let level = 0; level < errors.length; level++) {
  let total = 0;
  for (let i = 0; i < errors[level].length; i++) total += errors[level][i];
  means.push(total / errors[level].length);
  console.log(levelSteps[level] + ' steps: mean error ' + means[level].toFixed(6));
}
for (let level = 1; level < means.length; level++) {
  console.log('halving dt divided the error by ' + (means[level - 1] / means[level]).toFixed(3));
}
```

---

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

[Previous task](https://gpu.rocks/learn/ode-integrators-62f4a3ff/4.md) · [Next task](https://gpu.rocks/learn/ode-integrators-62f4a3ff/6.md)