# Four Slopes, Weighted

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

Classical RK4 samples the slope four times per step — once at the start, twice
at the middle (the second one using the first one's estimate), once at the end — and
combines them `(k₁ + 2·k₂ + 2·k₃ + k₄) / 6`. Four evaluations,
**fourth** order: halving `dt` divides the error by
*sixteen*.

Which is why this task's step sizes look absurd. The coarsest level crosses
`t = 0…5` in **five** steps — a full radian each — and it is
still more accurate than 400 Euler steps. Run RK4 at Euler's finest setting instead and
its error would be around 10⁻⁹, roughly a thousand times below what float32 can resolve
next to a value of order 1; the study would measure rounding, not convergence. Three
levels is what fits between "not yet asymptotic" and "already noise" — take one more
halving and the ratio stops being a clean 16.

## Goal

**Goal:** write the RK4 step and watch the ratio jump to ≈16.

## Requirements

- Four slope pairs: `k1` at the start, `k2` and `k3` at the midpoint, `k4` at the end
- `k2` and `k3` step out by `dt/2`; `k4` by the full `dt`
- Each probe builds on the *previous* one — `k3` from `k2`, `k4` from `k3`
- Combine with weights `1, 2, 2, 1` over `6`

## Hint 1 — the derivative of the state

The state is `(x, v)`, so every slope is a pair: the derivative
of `x` is `v`, and the derivative of `v` is
`−x`. That gives `k1x = v` and `k1v = −x`, and
every later probe is the same rule evaluated at a shifted state.

## Hint 2 — the middle two

`k2` is the slope half a step along `k1`;
`k3` is the slope half a step along `k2`:

```js
const k2x = v + k1v * half;
const k2v = -(x + k1x * half);
const k3x = v + k2v * half;
const k3v = -(x + k2x * half);
```

Building `k3` from `k1` again is the classic slip, and it
costs you two whole orders.

## Hint 3 — the combination

```js
x = x + (dt / 6) * (k1x + 2 * k2x + 2 * k3x + k4x);
v = v + (dt / 6) * (k1v + 2 * k2v + 2 * k3v + k4v);
```

The middle pair carries twice the weight, and the six is what the weights sum
to.

## Same idea elsewhere

RK4 trades memory traffic for arithmetic: four evaluations of the derivative per
step, all of them on values already sitting in registers. That is exactly the bargain GPUs
reward — arithmetic is nearly free, and the force evaluation in a real simulation is
usually memory-bound, so fewer, bigger, better steps often beat more cheap ones. It is the
same arithmetic-intensity argument that drives kernel fusion in CUDA, WebGPU and Metal.

## Starter code

```js
// Only THREE levels here, and they are brutally coarse: 5, 10, 20 steps
// to cross t = 0…5. Fourth order needs big steps to stay visible above float32.
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++) {
    const half = dt / 2;
    // TODO: four slope pairs — k1 at the start, k2 and k3 at the midpoint
    // (each from the previous one), k4 at the end — then combine them
    // with weights 1, 2, 2, 1 over 6.
    const k1x = v;
    const k1v = -x;
    x = x + k1x * dt;
    v = v + k1v * 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, 3],
  constants: { tEnd: 5 },
  loopMaxIterations: 20,
});

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/4

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