RK4 in JavaScript: A 30-Line Integrator

1. Why Euler Isn't Enough

The simplest way to integrate dy/dt = f(t, y) numerically is explicit (forward) Euler: take the current slope and step along it.

y += f(t, y) * dt;
t += dt;

It's one line, but its local error is O(dt²) per step — it only uses the derivative at the start of the interval, so any curvature in the true trajectory during that step is missed entirely. For a chaotic system like the Lorenz equations, that error doesn't just make the trajectory slightly wrong — because nearby trajectories diverge exponentially, a small integration error quickly becomes a completely different orbit. At larger step sizes Euler can also be unstable, with energy or amplitude drifting until the simulation visibly explodes.

2. The Four RK4 Stages

4th-order Runge-Kutta (RK4) samples the derivative four times per step — at the start, at two estimates of the midpoint, and at an estimate of the endpoint — then blends them with weights 1 : 2 : 2 : 1:

Intuitively: k1 is the slope you'd get from plain Euler. k2 and k3 refine that guess by evaluating the slope at the midpoint of the interval, using the previous stage's estimate to get there. k4 evaluates the slope at the far end, using the midpoint estimate to project forward. Averaging all four — with the midpoint estimates counted twice — cancels out the leading error terms, giving local error O(dt⁵) per step (global error O(dt⁴)) for only 4× the cost of one Euler step.

3. The 30-Line Integrator

The trick to a reusable RK4 implementation is to keep the state as a plain array and write two tiny vector helpers (add, scale) instead of hard-coding x, y, z. The stepper then works unchanged whether y has 2 components (a pendulum) or 3 (Lorenz) or 40 (a chain of coupled oscillators).

// --- vector helpers -------------------------------------------------
function add(a, b) {
  return a.map((v, i) => v + b[i]);
}
function scale(a, s) {
  return a.map((v) => v * s);
}

// --- one RK4 step -----------------------------------------------------
// f(t, y) -> array of derivatives, same shape as y
function rk4Step(f, t, y, dt) {
  const k1 = f(t,              y);
  const k2 = f(t + dt / 2, add(y, scale(k1, dt / 2)));
  const k3 = f(t + dt / 2, add(y, scale(k2, dt / 2)));
  const k4 = f(t + dt,     add(y, scale(k3, dt)));

  const sum = k1.map(
    (_, i) => k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]
  );
  return add(y, scale(sum, dt / 6));
}

// --- drive the simulation ---------------------------------------------
function integrate(f, y0, dt, steps) {
  let t = 0, y = y0, trail = [y0];
  for (let i = 0; i < steps; i++) {
    y = rk4Step(f, t, y, dt);
    t += dt;
    trail.push(y);
  }
  return trail;
}

Count it up: two 3-line vector helpers, an 11-line rk4Step, and an 8-line driver loop — right around 30 lines for a fully generic ODE integrator, with no dependencies.

4. Testing on the Lorenz System

To integrate the Lorenz system, all that's needed is a derivative function that matches the f(t, y) signature above — rk4Step itself doesn't change:

function lorenz(t, [x, y, z]) {
  const sigma = 10, rho = 28, beta = 8 / 3;
  return [
    sigma * (y - x),
    x * (rho - z) - y,
    x * y - beta * z,
  ];
}

const trail = integrate(lorenz, [0.1, 0, 0], 0.01, 20000);
// trail[i] = [x, y, z] at step i — feed straight into a WebGL line strip

Swap lorenz for a pendulum derivative ([θ, ω] → [ω, -g/L·sin θ]) or a Rössler system and the same rk4Step/integrate pair keeps working — that's the whole point of keeping state as a generic array.

5. Choosing a Step Size

RK4's local error scales as dt⁵, so halving dt cuts the per-step error to about 1/32 of its previous value — but doubles the number of steps needed to cover the same time span. In practice:

• dt = 0.005–0.01 works well for the Lorenz system at the classic parameters (σ=10, ρ=28, β=8/3).
• Check conservation, not just looks: for energy-conserving systems (an undamped pendulum), plot total energy over time — if it drifts noticeably, shrink dt.
• Fixed vs adaptive: this tutorial uses a fixed step. For production-grade accuracy with fewer wasted steps, an adaptive method (e.g. Dormand-Prince, RK45) adjusts dt automatically based on estimated local error — overkill for most interactive WebGL demos, but worth knowing about.

6. Common Pitfalls