Runge-Kutta Methods

Every time a video game drops a ball, a spacecraft planner traces a trajectory, or a weather model ticks forward an hour, a computer must answer the same question: given where things are now and the rule for how they change, where will they be a moment later?

That rule is a differential equation — it says "the velocity changes at this rate" or "the temperature cools at this rate." The trouble is that real rates depend on position, and position changes as you move. Follow the slope at the start and you overshoot; the path curves away from you.

The naive fix is Euler's method: take tiny steps, update position using the slope at the start of each step. Small enough steps and it works — but shrinking the step size costs proportionally more computation. You need to be smarter.

Runge-Kutta methods are that smarter idea. Instead of sampling the slope once per step, sample it several times — at the start, the midpoint, the end — then blend those samples into a single high-quality estimate. The classic fourth-order method (RK4), published by mathematicians Carl Runge and Wilhelm Kutta around 1900, gets the answer right to the fourth power of the step size. Halve the step, and the error drops by a factor of 16, not 2. It is the undisputed workhorse of numerical simulation.

The demo below integrates a simple circular orbit using both Euler's method and RK4 with the same step size. The planet starts at (1, 0) with velocity (0, 1), pulled toward the origin by gravity — the exact orbit is a perfect circle.

<div class="controls">
  <label>Step size h = <span id="h-val">0.3</span>
    <input type="range" id="h-slider" min="1" max="20" value="6" step="1">
  </label>
  <button id="run-btn" type="button">Run one orbit</button>
  <button id="reset-btn" type="button" class="ghost">Reset</button>
</div>
<canvas id="cvs" width="440" height="300"></canvas>
<div class="legend">
  <span class="dot euler"></span> Euler &nbsp;
  <span class="dot rk4"></span> RK4 &nbsp;
  <span class="dot exact"></span> Exact
</div>
<div id="errors" class="errors"></div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; background: #fff; }
.controls { display: flex; align-items: center; gap: .8rem; flex-wrap: wrap; margin-bottom: .5rem; font-size: .9rem; }
label { display: flex; align-items: center; gap: .4rem; }
input[type=range] { width: 120px; }
button { font: 600 13px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 7px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
canvas { display: block; border: 1px solid #dde3ea; border-radius: 8px; background: #f7f9fb; }
.legend { font-size: .82rem; margin: .4rem 0; display: flex; gap: 1rem; }
.dot { display: inline-block; width: 10px; height: 10px; border-radius: 50%; vertical-align: middle; }
.dot.euler { background: #e63946; }
.dot.rk4   { background: #2a9d8f; }
.dot.exact { background: #457b9d; }
.errors { font-size: .85rem; color: #333; min-height: 2.4em; }
.errors span { margin-right: 1rem; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Circular orbit: x'' = -x/r^3, y'' = -y/r^3  (GM=1, so circular at r=1 with v=1).
// State vector: [x, y, vx, vy]
// Derivatives: [vx, vy, -x/r^3, -y/r^3]

function deriv(s) {
  const r3 = Math.pow(s[0]*s[0] + s[1]*s[1], 1.5);
  return [s[2], s[3], -s[0]/r3, -s[1]/r3];
}

function eulerStep(s, h) {
  const d = deriv(s);
  return s.map((v, i) => v + h * d[i]);
}

function rk4Step(s, h) {
  const k1 = deriv(s);
  const s2 = s.map((v, i) => v + 0.5 * h * k1[i]);
  const k2 = deriv(s2);
  const s3 = s.map((v, i) => v + 0.5 * h * k2[i]);
  const k3 = deriv(s3);
  const s4 = s.map((v, i) => v + h * k3[i]);
  const k4 = deriv(s4);
  return s.map((v, i) => v + (h / 6) * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]));
}

function integrateOrbit(stepFn, h) {
  const T = 2 * Math.PI; // one full orbit
  const steps = Math.ceil(T / h);
  const pts = [];
  let s = [1, 0, 0, 1]; // start: position (1,0), velocity (0,1)
  pts.push([s[0], s[1]]);
  for (let i = 0; i < steps; i++) {
    s = stepFn(s, h);
    pts.push([s[0], s[1]]);
  }
  return { pts, final: s };
}

function exactPos(t) {
  return [Math.cos(t), Math.sin(t)];
}

const cvs = document.getElementById('cvs');
const ctx = cvs.getContext('2d');
const W = cvs.width, H = cvs.height;
const cx = W / 2, cy = H / 2;
const scale = 110; // pixels per unit

function toCanvas(x, y) {
  return [cx + x * scale, cy - y * scale];
}

function drawPath(pts, color, lineWidth) {
  ctx.beginPath();
  ctx.strokeStyle = color;
  ctx.lineWidth = lineWidth;
  const [x0, y0] = toCanvas(pts[0][0], pts[0][1]);
  ctx.moveTo(x0, y0);
  for (let i = 1; i < pts.length; i++) {
    const [xi, yi] = toCanvas(pts[i][0], pts[i][1]);
    ctx.lineTo(xi, yi);
  }
  ctx.stroke();
}

function drawDot(x, y, color, r) {
  const [cx2, cy2] = toCanvas(x, y);
  ctx.beginPath();
  ctx.arc(cx2, cy2, r, 0, 2 * Math.PI);
  ctx.fillStyle = color;
  ctx.fill();
}

function drawExactCircle() {
  const [ex, ey] = toCanvas(0, 0);
  ctx.beginPath();
  ctx.arc(ex, ey, scale, 0, 2 * Math.PI);
  ctx.strokeStyle = '#457b9d';
  ctx.lineWidth = 1.5;
  ctx.setLineDash([5, 4]);
  ctx.stroke();
  ctx.setLineDash([]);
}

function getH() {
  const slider = document.getElementById('h-slider');
  // map slider 1..20 -> h 0.05..1.5
  return 0.05 + (parseInt(slider.value) - 1) * (1.45 / 19);
}

function run() {
  const h = getH();
  const euler = integrateOrbit(eulerStep, h);
  const rk = integrateOrbit(rk4Step, h);
  const T = 2 * Math.PI;
  const [ex, ey] = exactPos(T); // = (1,0) after full period

  // Errors: distance from exact position after one orbit
  const eulerErr = Math.hypot(euler.final[0] - ex, euler.final[1] - ey);
  const rkErr    = Math.hypot(rk.final[0] - ex,    rk.final[1] - ey);

  ctx.clearRect(0, 0, W, H);
  drawExactCircle();
  drawPath(euler.pts, '#e63946', 1.8);
  drawPath(rk.pts, '#2a9d8f', 1.8);

  // Draw planet dots at final position
  drawDot(euler.final[0], euler.final[1], '#e63946', 5);
  drawDot(rk.final[0], rk.final[1], '#2a9d8f', 5);
  // Draw sun at origin
  drawDot(0, 0, '#f4a261', 7);

  // Log for node verification
  console.log(fmt('h=%h% euler_err=%eulerErr% rk4_err=%rkErr%', { h: h.toFixed(3), eulerErr: eulerErr.toFixed(6), rkErr: rkErr.toFixed(6) }));

  const errDiv = document.getElementById('errors');
  errDiv.innerHTML =
    fmt('<span><b style="color:#e63946">Euler</b> error: %eulerErr%</span><span><b style="color:#2a9d8f">RK4</b> error: %rkErr%</span>', { eulerErr: eulerErr.toFixed(4), rkErr: rkErr.toFixed(4) });
}

const slider = document.getElementById('h-slider');
const hVal = document.getElementById('h-val');
slider.addEventListener('input', () => {
  hVal.textContent = getH().toFixed(2);
});
document.getElementById('run-btn').addEventListener('click', run);
document.getElementById('reset-btn').addEventListener('click', () => {
  ctx.clearRect(0, 0, W, H);
  document.getElementById('errors').innerHTML = '';
});

// Initial draw
hVal.textContent = getH().toFixed(2);
run();

<div class="controls">
  <label>Step size h = <span id="h-val">0.3</span>
    <input type="range" id="h-slider" min="1" max="20" value="6" step="1">
  </label>
  <button id="run-btn" type="button">Run one orbit</button>
  <button id="reset-btn" type="button" class="ghost">Reset</button>
</div>
<canvas id="cvs" width="440" height="300"></canvas>
<div class="legend">
  <span class="dot euler"></span> Euler &nbsp;
  <span class="dot rk4"></span> RK4 &nbsp;
  <span class="dot exact"></span> Exact
</div>
<div id="errors" class="errors"></div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; background: #fff; }
.controls { display: flex; align-items: center; gap: .8rem; flex-wrap: wrap; margin-bottom: .5rem; font-size: .9rem; }
label { display: flex; align-items: center; gap: .4rem; }
input[type=range] { width: 120px; }
button { font: 600 13px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 7px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
canvas { display: block; border: 1px solid #dde3ea; border-radius: 8px; background: #f7f9fb; }
.legend { font-size: .82rem; margin: .4rem 0; display: flex; gap: 1rem; }
.dot { display: inline-block; width: 10px; height: 10px; border-radius: 50%; vertical-align: middle; }
.dot.euler { background: #e63946; }
.dot.rk4   { background: #2a9d8f; }
.dot.exact { background: #457b9d; }
.errors { font-size: .85rem; color: #333; min-height: 2.4em; }
.errors span { margin-right: 1rem; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Circular orbit: x'' = -x/r^3, y'' = -y/r^3  (GM=1, so circular at r=1 with v=1).
// State vector: [x, y, vx, vy]
// Derivatives: [vx, vy, -x/r^3, -y/r^3]

function deriv(s) {
  const r3 = Math.pow(s[0]*s[0] + s[1]*s[1], 1.5);
  return [s[2], s[3], -s[0]/r3, -s[1]/r3];
}

function eulerStep(s, h) {
  const d = deriv(s);
  return s.map((v, i) => v + h * d[i]);
}

function rk4Step(s, h) {
  const k1 = deriv(s);
  const s2 = s.map((v, i) => v + 0.5 * h * k1[i]);
  const k2 = deriv(s2);
  const s3 = s.map((v, i) => v + 0.5 * h * k2[i]);
  const k3 = deriv(s3);
  const s4 = s.map((v, i) => v + h * k3[i]);
  const k4 = deriv(s4);
  return s.map((v, i) => v + (h / 6) * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]));
}

function integrateOrbit(stepFn, h) {
  const T = 2 * Math.PI; // one full orbit
  const steps = Math.ceil(T / h);
  const pts = [];
  let s = [1, 0, 0, 1]; // start: position (1,0), velocity (0,1)
  pts.push([s[0], s[1]]);
  for (let i = 0; i < steps; i++) {
    s = stepFn(s, h);
    pts.push([s[0], s[1]]);
  }
  return { pts, final: s };
}

function exactPos(t) {
  return [Math.cos(t), Math.sin(t)];
}

const cvs = document.getElementById('cvs');
const ctx = cvs.getContext('2d');
const W = cvs.width, H = cvs.height;
const cx = W / 2, cy = H / 2;
const scale = 110; // pixels per unit

function toCanvas(x, y) {
  return [cx + x * scale, cy - y * scale];
}

function drawPath(pts, color, lineWidth) {
  ctx.beginPath();
  ctx.strokeStyle = color;
  ctx.lineWidth = lineWidth;
  const [x0, y0] = toCanvas(pts[0][0], pts[0][1]);
  ctx.moveTo(x0, y0);
  for (let i = 1; i < pts.length; i++) {
    const [xi, yi] = toCanvas(pts[i][0], pts[i][1]);
    ctx.lineTo(xi, yi);
  }
  ctx.stroke();
}

function drawDot(x, y, color, r) {
  const [cx2, cy2] = toCanvas(x, y);
  ctx.beginPath();
  ctx.arc(cx2, cy2, r, 0, 2 * Math.PI);
  ctx.fillStyle = color;
  ctx.fill();
}

function drawExactCircle() {
  const [ex, ey] = toCanvas(0, 0);
  ctx.beginPath();
  ctx.arc(ex, ey, scale, 0, 2 * Math.PI);
  ctx.strokeStyle = '#457b9d';
  ctx.lineWidth = 1.5;
  ctx.setLineDash([5, 4]);
  ctx.stroke();
  ctx.setLineDash([]);
}

function getH() {
  const slider = document.getElementById('h-slider');
  // map slider 1..20 -> h 0.05..1.5
  return 0.05 + (parseInt(slider.value) - 1) * (1.45 / 19);
}

function run() {
  const h = getH();
  const euler = integrateOrbit(eulerStep, h);
  const rk = integrateOrbit(rk4Step, h);
  const T = 2 * Math.PI;
  const [ex, ey] = exactPos(T); // = (1,0) after full period

  // Errors: distance from exact position after one orbit
  const eulerErr = Math.hypot(euler.final[0] - ex, euler.final[1] - ey);
  const rkErr    = Math.hypot(rk.final[0] - ex,    rk.final[1] - ey);

  ctx.clearRect(0, 0, W, H);
  drawExactCircle();
  drawPath(euler.pts, '#e63946', 1.8);
  drawPath(rk.pts, '#2a9d8f', 1.8);

  // Draw planet dots at final position
  drawDot(euler.final[0], euler.final[1], '#e63946', 5);
  drawDot(rk.final[0], rk.final[1], '#2a9d8f', 5);
  // Draw sun at origin
  drawDot(0, 0, '#f4a261', 7);

  // Log for node verification
  console.log(fmt('h=%h% euler_err=%eulerErr% rk4_err=%rkErr%', { h: h.toFixed(3), eulerErr: eulerErr.toFixed(6), rkErr: rkErr.toFixed(6) }));

  const errDiv = document.getElementById('errors');
  errDiv.innerHTML =
    fmt('<span><b style="color:#e63946">Euler</b> error: %eulerErr%</span><span><b style="color:#2a9d8f">RK4</b> error: %rkErr%</span>', { eulerErr: eulerErr.toFixed(4), rkErr: rkErr.toFixed(4) });
}

const slider = document.getElementById('h-slider');
const hVal = document.getElementById('h-val');
slider.addEventListener('input', () => {
  hVal.textContent = getH().toFixed(2);
});
document.getElementById('run-btn').addEventListener('click', run);
document.getElementById('reset-btn').addEventListener('click', () => {
  ctx.clearRect(0, 0, W, H);
  document.getElementById('errors').innerHTML = '';
});

// Initial draw
hVal.textContent = getH().toFixed(2);
run();

Drag the step size slider to make steps larger. With a large step, Euler's planet spirals outward — energy leaks in with every step — while RK4 holds its orbit. Shrink the step and both methods improve, but RK4 converges far faster. The error printed below the canvas is the distance between each planet and the exact position after one full orbit.

Runge-Kutta is not a mystery or an open problem — it is a solved algorithm with proven, quantified accuracy. Here is what the analysis says:

• Euler's method is order 1: the local error per step is proportional to $h^{2}$, so the accumulated error over a fixed time interval is $O(h)$. Halve the step, halve the error.
• RK4 is order 4: local error is $O(h^{5})$, accumulated error is $O(h^{4})$. Halve the step, drop the error by a factor of 16. Same interval, same accuracy, far fewer steps.
• The RK4 recipe at each step: evaluate the slope $k_{1}$ at the start, $k_{2}$ at the midpoint using $k_{1}$ to advance, $k_{3}$ at the midpoint using $k_{2}$, and $k_{4}$ at the end using $k_{3}$. Blend them as ($k_{1}$ + $2k_{2}$ + $2k_{3}$ + $k_{4}$) / 6. That weighted average is the optimal estimate for a fourth-order method.
• Status: solved. RK4 was fully analyzed by the early twentieth century. It is not NP-complete, not undecidable, not an open problem. You can compute the error bound analytically for any smooth ODE.
• Adaptive variants (Dormand–Prince, Runge–Kutta–Fehlberg) run two methods of different orders simultaneously, estimate the local error, and automatically shrink or grow the step size to hit a target tolerance — this is what scipy.integrate.solve_ivp and MATLAB's ode45 implement under the hood.

The only "hard" part of numerical integration is stiffness: some ODEs have components that change on wildly different time scales, forcing very small steps to stay stable. Stiff solvers (implicit methods such as Radau or BDF) handle these at the cost of solving a small system of equations per step, but RK4 remains the right tool for the vast majority of smooth, non-stiff problems.

Wherever a system's future depends on its current state through a rate law, Runge-Kutta is in the room:

• Orbital mechanics and spaceflight: trajectory planning for satellites, the Moon, and interplanetary probes all integrate Newton's law of gravitation with RK4 or adaptive variants. The N-body simulation problem is exactly this, scaled up.
• Game and film physics engines: rigid-body dynamics, cloth, fluids, and particle systems in real-time engines use RK4 (or symplectic variants that conserve energy) to update positions and velocities each frame.
• Climate and weather models: differential equations for temperature, pressure, and humidity are stepped forward in time — typically with semi-implicit schemes that combine RK4's accuracy with stiff-solver stability.
• Electrical circuit simulation: SPICE-type simulators integrate the ODEs of circuit nodes; Dormand–Prince is the default integrator in many such tools.
• Pharmacokinetics: drug concentration in blood compartments follows first-order ODEs, and RK4 lets researchers model dosing regimens without a closed-form solution.
• Machine learning: the ResNet architecture can be interpreted as an Euler discretization of an ODE, and Neural ODEs (Chen et al., 2018) explicitly replace discrete layers with a continuous ODE solved by an adaptive RK integrator — bridging deep learning and numerical analysis.

Runge-Kutta methods are a masterclass in turning a naive idea into a precise one. Euler's method says "follow the slope from where you stand." RK4 says "sample the slope at four carefully chosen points, then follow their weighted average." The cost is three extra function evaluations per step. The reward is an error that shrinks sixteen times faster as you refine the step.

Carl Runge and Wilhelm Kutta published the core idea around 1900 — long before computers existed — and it has survived every decade of numerical analysis since. Modern adaptive solvers have refined the packaging, but the heart of the algorithm is unchanged: measure the landscape at several points before you step.

The next time your favorite game drops a physics object, a satellite stays in orbit, or a simulation unfolds without numerical drift, a variation of this idea is probably running underneath. Four slope samples, a weighted blend, and a century of trust — that is all it takes to simulate the universe on a digital machine.