Introduction

Every time a physicist simulates a satellite orbit, a biologist models a chemical reaction, or an engineer tests a control loop, they are solving a differential equation: given how fast a quantity changes now, trace where it goes. No closed formula exists for most real problems, so we march forward in small time steps, computing an approximate value at each one.

The naive approach picks a fixed step size: the same tiny increment from start to finish. That is wasteful. Most trajectories have long, gentle stretches where almost nothing happens — and short, violent bursts where the state changes rapidly. A fixed step that is small enough to survive the burst spends the vast majority of its budget plodding through the calm.

Adaptive step-size control solves both problems at once. The solver estimates its own local error at each step — the difference between a cheap approximation and a slightly better one. If the error is too big, it shrinks the step and retries. If the error is tiny, it stretches the step and moves faster. The result is a solver that automatically clusters its effort at the sharp turns and coasts through the smooth stretches.

This is the core idea behind workhorses like Runge-Kutta 4(5) (Dormand–Prince), the algorithm that powers SciPy's solve_ivp, MATLAB's ode45, and countless scientific applications. The technique is decades old and fully understood — it is not an open problem but a solved piece of numerical craft worth knowing.

Try It

The demo below integrates y′=2sin⁡(2t)+5 e−50(t−1)2y' = 2\sin(2t) + 5\,e^{-50(t-1)^2} on [0,6][0, 6]: a smooth sine wave carrying a sharp Gaussian spike near t=1t = 1. A fixed-step Euler solver and a simple adaptive RK4 solver run side by side — both chasing the same curve.

<!-- Demo layout for adaptive step-size integration article -->
<p class="hint">Both solvers integrate <i>y</i>&#x27; = 2&thinsp;sin(2<i>t</i>) + 5&thinsp;exp(&minus;50(<i>t</i>&minus;1)&#xB2;) — a smooth wave with a <b>sharp Gaussian spike</b> near <i>t</i>&thinsp;=&thinsp;1. <b>Gray dots</b> = fixed-step Euler (40 steps, evenly spaced). <b>Red dots</b> = adaptive RK4 — watch them cluster at the spike and spread out elsewhere. Drag the slider to change the tolerance.</p>
<div class="controls">
  <label>Tolerance: <input id="tol" type="range" min="1" max="5" step="1" value="3"></label>
  <span id="tol-display" class="tol-val"></span>
</div>
<canvas id="cv" width="560" height="260"></canvas>
<div class="legend">
  <span class="dot fixed"></span> Fixed-step Euler (40 steps)
  &nbsp;&nbsp;
  <span class="dot adaptive"></span> Adaptive RK4
</div>
<div class="stats" id="stats"></div>
<div class="btns">
  <button id="run" type="button">Run solvers</button>
  <button id="reset" type="button" class="ghost">Reset</button>
</div>

/* Styles for the adaptive step-size demo */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; }
.hint { font-size: .9rem; color: #444; margin: 0 0 .6rem; line-height: 1.45; }
.controls { display: flex; align-items: center; gap: .7rem; margin-bottom: .5rem; font-size: .9rem; }
input[type=range] { width: 120px; }
.tol-val { font-weight: 700; font-size: .9rem; color: #1d3557; min-width: 4rem; }
canvas { display: block; width: 100%; max-width: 560px; border: 1px solid #cdd9e3; border-radius: 8px; background: #f8fafc; }
.legend { display: flex; align-items: center; gap: .5rem; font-size: .85rem; margin: .4rem 0; flex-wrap: wrap; }
.dot { display: inline-block; width: 10px; height: 10px; border-radius: 50%; }
.dot.fixed { background: #adb1b8; }
.dot.adaptive { background: #e63946; }
.stats { font-size: .9rem; font-weight: 600; min-height: 1.4em; color: #1d3557; margin: .2rem 0; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; margin-top: .4rem; }
button { font: 600 14px system-ui, sans-serif; padding: .45rem .9rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }

/* Adaptive step-size integration demo — fixed Euler vs adaptive RK4 */
const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));

/* ODE right-hand side (autonomous in t only, not y) */
/* y' = 2*sin(2*t) + 5*exp(-50*(t-1)^2), y(0) = 0 */
/* smooth sine wave plus a sharp Gaussian spike near t = 1 */
function f(t) { return 2 * Math.sin(2 * t) + 5 * Math.exp(-50 * (t - 1) * (t - 1)); }

const T_END = 6;
const Y0 = 0; /* initial condition */

/* tolerance levels mapped to slider positions 1 through 5 */
const TOL_LEVELS = [1e-1, 5e-2, 1e-2, 1e-3, 1e-4];
const TOL_LABELS = ['1e-1', '5e-2', '1e-2', '1e-3', '1e-4'];

/* Fixed-step Euler integrator — same step size everywhere */
function fixedEuler(n) {
  const h = T_END / n;
  const pts = [{ t: 0, y: Y0 }];
  let t = 0, y = Y0;
  for (let i = 0; i < n; i++) { y += h * f(t); t += h; pts.push({ t, y }); }
  return pts;
}

/* Adaptive RK4 with embedded error estimate and step rescaling */
function adaptiveRK(tol) {
  const pts = [{ t: 0, y: Y0 }];
  let t = 0, y = Y0, h = 0.3;
  const hMin = 1e-7, hMax = 1.0, safety = 0.9, order = 4;
  let steps = 0;
  while (t < T_END && steps < 50000) {
    if (t + h > T_END) h = T_END - t;
    /* fourth-order Runge-Kutta estimate */
    const k1 = f(t);
    const k2 = f(t + h / 2);
    const k3 = f(t + h / 2);
    const k4 = f(t + h);
    const y4 = y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
    /* fifth-order correction used only for error estimation */
    const k5 = f(t + h);
    const y5 = y + h / 90 * (7 * k1 + 32 * k3 + 12 * k4 + 32 * k5 + 7 * k1);
    const err = Math.abs(y5 - y4);
    const scale = tol + tol * Math.abs(y4);
    if (err <= scale || h <= hMin) {
      t += h; y = y4; pts.push({ t, y }); /* accept step and record point */
    }
    /* rescale h based on error-to-tolerance ratio */
    const factor = err > 0 ? safety * Math.pow(scale / err, 1 / (order + 1)) : 2;
    h = Math.min(hMax, Math.max(hMin, h * factor));
    steps++;
  }
  return pts;
}

const canvas = document.getElementById('cv');
const ctx = canvas.getContext('2d');
const W = canvas.width, H = canvas.height;
const PAD = { top: 18, right: 18, bottom: 28, left: 42 };
const Y_MIN = -0.5, Y_MAX = 2.5;

function toPixel(t, y) {
  const x = PAD.left + (t / T_END) * (W - PAD.left - PAD.right);
  const yNorm = (y - Y_MIN) / (Y_MAX - Y_MIN);
  const py = PAD.top + (1 - yNorm) * (H - PAD.top - PAD.bottom);
  return { x, py };
}

function drawAxes() {
  /* horizontal grid lines and vertical time ticks */
  ctx.save();
  ctx.strokeStyle = '#cdd9e3'; ctx.lineWidth = 1;
  for (let v = 0; v <= 2; v += 0.5) {
    const { py } = toPixel(0, v);
    ctx.beginPath(); ctx.moveTo(PAD.left, py); ctx.lineTo(W - PAD.right, py); ctx.stroke();
    ctx.fillStyle = '#666'; ctx.font = '10px system-ui'; ctx.textAlign = 'right';
    ctx.fillText(v.toFixed(1), PAD.left - 4, py + 3);
  }
  for (let tv = 0; tv <= T_END + 0.01; tv += 1) {
    const { x } = toPixel(tv, 0);
    ctx.beginPath(); ctx.moveTo(x, PAD.top); ctx.lineTo(x, H - PAD.bottom); ctx.stroke();
    ctx.fillStyle = '#666'; ctx.textAlign = 'center';
    ctx.fillText(tv.toFixed(0), x, H - PAD.bottom + 12);
  }
  ctx.restore();
}

function drawExact() {
  /* draw reference trajectory with fine fixed steps */
  ctx.save(); ctx.strokeStyle = '#1d355740'; ctx.lineWidth = 1.5;
  ctx.beginPath();
  let y = Y0;
  for (let i = 0; i <= 600; i++) {
    const t = (i / 600) * T_END;
    if (i > 0) { const dt = T_END / 600; y += dt * f(t - dt); }
    const { x, py } = toPixel(t, y);
    i === 0 ? ctx.moveTo(x, py) : ctx.lineTo(x, py);
  }
  ctx.stroke(); ctx.restore();
}

function drawLine(pts, color) {
  ctx.save(); ctx.strokeStyle = color; ctx.lineWidth = 1.5;
  ctx.beginPath();
  pts.forEach(({ t, y }, i) => {
    const { x, py } = toPixel(t, y);
    i === 0 ? ctx.moveTo(x, py) : ctx.lineTo(x, py);
  });
  ctx.stroke(); ctx.restore();
}

function drawDots(pts, color, r) {
  /* draw each accepted solver step as a dot */
  ctx.save(); ctx.fillStyle = color;
  pts.forEach(({ t, y }) => {
    const { x, py } = toPixel(t, y);
    ctx.beginPath(); ctx.arc(x, py, r, 0, 2 * Math.PI); ctx.fill();
  });
  ctx.restore();
}

const N_FIXED = 40; /* fixed number of Euler steps for comparison */

function run() {
  const tolIdx = parseInt(document.getElementById('tol').value, 10) - 1;
  const tol = TOL_LEVELS[tolIdx];
  const fixedPts = fixedEuler(N_FIXED);
  const adaptPts = adaptiveRK(tol);
  ctx.clearRect(0, 0, W, H);
  drawAxes();
  drawExact();
  drawLine(fixedPts, '#adb1b8');
  drawLine(adaptPts, '#e63946');
  drawDots(fixedPts, '#8a8f99', 3);
  drawDots(adaptPts, '#e63946', 3.5);
  document.getElementById('stats').textContent =
    i18nFmt('Fixed: %fixed% points · Adaptive: %adaptive% points · tol = %tol%', { fixed: N_FIXED + 1, adaptive: adaptPts.length, tol: TOL_LABELS[tolIdx] });
}

function updateTolDisplay() {
  const idx = parseInt(document.getElementById('tol').value, 10) - 1;
  document.getElementById('tol-display').textContent = 'tol = ' + TOL_LABELS[idx];
}

document.getElementById('tol').addEventListener('input', updateTolDisplay);
document.getElementById('run').addEventListener('click', run);
document.getElementById('reset').addEventListener('click', () => {
  document.getElementById('tol').value = 3;
  updateTolDisplay();
  ctx.clearRect(0, 0, W, H);
  drawAxes();
  drawExact();
  document.getElementById('stats').textContent = '';
});

/* initialize slider label and draw on page load */
updateTolDisplay();
run();
Result
Console

Watch how the adaptive solver's dots cluster around the spike near t=1t = 1 and then spread far apart on the smooth sinusoidal tail. The fixed solver uses the same number of steps everywhere — paying the same cost for the easy stretches as for the hard ones. The step-count gap shows the computational savings, and tightening the tolerance reveals how the adaptive solver automatically re-clusters.

The Real Complexity

The machinery inside an adaptive solver is surprisingly compact.

Two approximations per step. A solver like RK4(5) computes two estimates of the next state: one of order pp and one of order p+1p+1. Their difference δ\delta is the local truncation error estimate — the cost of taking this step at this size.

The rescaling rule. If the current step hh produces error δ\delta but the user's tolerance is ε\varepsilon, the optimal next step is approximately:

hnew=h⋅(εδ)1/(p+1)h_{\text{new}} = h \cdot \left(\frac{\varepsilon}{\delta}\right)^{1/(p+1)}

This formula follows from the fact that the local error scales as hp+1h^{p+1}. Shrink hh by a factor of two and the error drops by 2p+12^{p+1} — roughly 32× for a fifth-order method.

The accept/reject loop. If δ>ε\delta > \varepsilon, the step is rejected and retried with hnewh_{\text{new}}. If δ≤ε\delta \leq \varepsilon, the step is accepted and the solver advances. A safety factor (typically 0.9) prevents oscillation around the boundary.

Stiffness is the hard case. Some ODEs — like chemical reaction networks with very different time scales — force the step size so small that even an adaptive explicit solver crawls. That is the domain of implicit methods and the separate field of stiff equation solvers, which invert a Jacobian at each step to buy larger steps at the cost of linear algebra.

The bottom line: adaptive step-size control is O(S)O(S) in the number of accepted steps SS, where SS is determined automatically by the solution's curvature and the user's tolerance. A factor-of-10 tighter tolerance roughly multiplies SS by 101/(p+1)10^{1/(p+1)} — for a fifth-order method, less than a factor of 2. Tight tolerances are cheaper than they look.

Where It Matters

Any simulation that follows a changing quantity through time benefits from step adaptivity:

  • Orbital mechanics: a satellite glides smoothly far from Earth but accelerates sharply during a close flyby. An adaptive solver naturally concentrates steps near the periapsis and coasts through the cruise phase, delivering accurate trajectories at a fraction of the cost of a fixed-step integrator.
  • Pharmacokinetics: drug concentration in the bloodstream drops fast after injection, then decays slowly. Adaptive solvers let pharmacologists fit models to clinical data without hand-tuning a step size for each drug.
  • Circuit simulation (SPICE): digital switching events create near-discontinuities in voltage and current. SPICE-family simulators use adaptive time-stepping to catch each edge accurately while flying through the quiet intervals.
  • Climate and weather models: atmospheric chemistry involves species with reaction rates spanning twelve orders of magnitude. Splitting the system and using adaptive solvers for the stiff chemical subsystem is standard practice.
  • Robotics and control: real-time simulators for legged robots must track contact events — sudden impacts — accurately. Adaptive solvers with event detection locate the exact moment of contact without sampling blindly.

The idea even migrated into numerical optimization: line-search methods and trust-region methods adaptively resize their "step" toward the minimum using exactly the same intuition — measure local curvature, adjust step size accordingly.

Conclusion

Adaptive step-size integration is a masterclass in computational economy. Instead of committing to a fixed budget per unit of time, the solver listens to the problem: it estimates how hard each local piece is and allocates effort accordingly. Sharp turns get many small steps; smooth stretches get a few large ones.

The rescaling formula hnew=h⋅(ε/δ)1/(p+1)h_{\text{new}} = h \cdot (\varepsilon / \delta)^{1/(p+1)} is a closed-loop controller — the algorithm controls its own accuracy in real time. Decades of software engineering have layered dense-output interpolation, event detection, and stiffness detection on top of this simple core, but the heart of it has not changed since the 1960s.

The lesson transfers broadly: wherever a computation has regions of wildly varying difficulty, the right strategy is almost always to measure difficulty locally and adapt. That idea appears in mesh refinement for numerical integration, in variable-rate sampling for signal processing, and in the step-size rules of gradient descent. The pattern recurs because it is simply optimal.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/adaptive-step-integration/Content licensed under CC BY-NC 4.0.