Adaptive Methods for ODEs

Steps an ODE solver with an adaptive step size that shrinks or grows to hold the local error near a tolerance.

Author

Apurva Nakade

Published

July 13, 2026

// VM is the shared utility library, loaded globally via _includes/head-scripts.html (js/**). This page uses VM.expressions.makeFunction2, VM.expressions.makeNumber, VM.plotting.paddedRange, VM.ui.renderTable, VM.numerical.rk4Solve, VM.ui.applyExampleParams.
VM = window.VM

methodOptions = ["Reference solution", "Adaptive RK45"]

// Matches the line/marker colors mainPlot uses for each trace, so the
// overlay legend's swatches double as a color key.
// vmTheme is what makes the colors below reactive. The chart's own chrome
// (background, axes, gridlines) repaints on a dark-mode toggle by itself,
// via the relayout patch in plotly-fullscreen-button.js -- but the *trace*
// colors are values this page bakes into the traces it builds, so nothing
// repaints them unless mainPlot re-runs. Depending on vmTheme is what
// re-runs it. See VM.plotting.onThemeChange for why it fires a couple of
// frames after the toggle rather than inline.
vmTheme = Generators.observe(notify => VM.plotting.onThemeChange(notify))

chartColors = VM.plotting.colors(vmTheme)

methodColors = new Map([
  ["Reference solution", chartColors.fn],
  ["Adaptive RK45", chartColors.alt]
])

// examples backs the "Try an example" dropdown near the top controls.
examples = [
  {
    title: "Fast attraction to cos(t)",
    params: {max: "2000", f: "-10 * (y - cos(t)) - sin(t)", y0: "0", b: "3"}
  },
  {
    title: "Stiff-looking decay",
    params: {max: "2000", f: "-15 * y", y0: "1", b: "1"}
  },
  {
    title: "Nonlinear Riccati-type equation",
    params: {max: "2000", f: "y^2 - t", y0: "0", b: "2"}
  },
  {
    title: "Oscillatory forcing with growth",
    params: {max: "2000", f: "y + sin(5 * t)", y0: "0", b: "3"}
  },
  {
    title: "Damped sinusoidal forcing",
    params: {max: "2000", f: "sin(t) - 0.5 * y", y0: "0", b: "10"}
  },
  {
    title: "Logistic growth",
    params: {max: "2000", f: "y * (1 - y)", y0: "0.1", b: "8"}
  }
]

viewof fText = {
  const params = new URLSearchParams(window.location.search)
  const fTextInput = Inputs.text({label: "y' = f(t, y)", value: params.get("f") ?? "-10 * (y - cos(t)) - sin(t)", placeholder: "Example: y - t^2 + 1, t - y, y * (1 - y)"})
  fTextInput.classList.add("ojs-fill")
  fTextInput.dataset.exampleField = "f"
  return fTextInput
}

// Plot commits the current function/initial-condition/endpoint fields —
// edits above don't take effect until this is clicked, so result/mainPlot
// don't recompute on every keystroke. The tolerance slider overlaid on the
// chart is unaffected and stays reactive.
viewof plotTrigger = {
  const plotButton = Inputs.button("Plot", {value: 0, reduce: v => v + 1})
  plotButton.classList.add("ojs-auto")
  return plotButton
}

viewof y0 = {
  const params = new URLSearchParams(window.location.search)
  const y0Input = Inputs.text({label: "Initial condition y(0)", value: params.get("y0") ?? "0", placeholder: "Example: 0.5, 1/2"})
  y0Input.dataset.exampleField = "y0"
  return y0Input
}

viewof bEnd = {
  const params = new URLSearchParams(window.location.search)
  const bEndInput = Inputs.text({label: "Right endpoint b", value: params.get("b") ?? "3", placeholder: "Example: 2, pi"})
  bEndInput.dataset.exampleField = "b"
  return bEndInput
}

// Always-live (not gated behind Plot) — it bounds the adaptive stepper
// rather than feeding result, so there's no reason to wait for a commit.
// Shares this row with the commit-gated fields now that a shorter label
// leaves room; see urlSyncControls below for how it stays in the URL.
viewof maxStepsInput = {
  const input = Inputs.text({value: initialMax, label: "Max steps", placeholder: "Example: 2000, 1000*2"})
  input.dataset.exampleField = "max"
  return input
}

viewof exampleSelect = {
  const options = [null, ...examples]
  const select = Inputs.select(options, {label: "Try an example", format: e => e ? e.title : "Pick an example…", value: null})
  select.classList.add("ojs-fill")
  return select
}

// Plain CSS selectors, not direct viewof references -- see
// js/ui/apply-example.js's comment for why.
exampleFieldSelectors = ({
  f: '[data-example-field="f"]',
  y0: '[data-example-field="y0"]',
  b: '[data-example-field="b"]',
  max: '[data-example-field="max"]'
})

// Applies the selected example's params to the fields above and clicks
// Plot, in place — no page reload, so no scroll jump either. A plain
// "change" listener attached once to the stable native <select> (rather
// than an OJS cell reactive on exampleSelect's value) — a <select> fires
// both "input" and "change" per selection, and a reactive cell would run
// once per event, with two overlapping calls racing each other.
applyExampleFromSelect = {
  const nativeSelect = (viewof exampleSelect).querySelector("select")
  if (!nativeSelect) return
  nativeSelect.addEventListener("change", () => {
    const example = (viewof exampleSelect).value
    if (!example) return
    VM.ui.applyExampleParams(exampleFieldSelectors, example.params, viewof plotTrigger)
  })
}

// Pressing Enter in the function/initial-condition/endpoint fields clicks
// Plot instead of submitting the field's own (invisible) form. Depends
// only on the stable views, so this wiring runs once and is never
// re-attached.
enterToPlot = {
  for (const view of [viewof fText, viewof y0, viewof bEnd]) {
    const input = view.querySelector("input")
    if (!input) continue
    input.addEventListener("keydown", (e) => {
      if (e.key !== "Enter") return
      e.preventDefault()
      viewof plotTrigger.querySelector("button").click()
    })
  }
}

// committed snapshots fText/y0/bEnd only when Plot is clicked, and writes
// them back into the URL so the address bar stays a shareable link. It
// reads `viewof fText`/`viewof y0`/`viewof bEnd` (the stable DOM views)
// rather than `fText`/`y0`/`bEnd` (the reactive values), so it does NOT
// rerun on every keystroke — only when `plotTrigger` changes. It still
// runs once on initial page load, so the first render uses the
// URL-provided (or default) values with no click.
committed = {
  plotTrigger
  const math = window.math
  const fVal = (viewof fText).value
  const y0Text = (viewof y0).value
  const bText = (viewof bEnd).value
  const y0Val = VM.expressions.makeNumber(math, y0Text)
  const bVal = VM.expressions.makeNumber(math, bText)

  const params = new URLSearchParams(window.location.search)
  params.set("f", fVal)
  params.set("y0", y0Text)
  params.set("b", bText)
  history.replaceState(null, "", `${window.location.pathname}?${params}`)

  return {fText: fVal, y0: y0Val, bEnd: bVal}
}

// max is always-live rather than commit-gated, so it gets its own sync cell
// instead of piggybacking on committed above (which only runs when Plot is
// clicked).
urlSyncControls = {
  const params = new URLSearchParams(window.location.search)
  params.set("max", maxStepsInput)
  history.replaceState(null, "", `${window.location.pathname}?${params}`)
  return true
}

// rk45Solve(f, t0, y0, tEnd, tol, maxSteps) advances y' = f(t, y), y(t0) = y0
// from t0 to tEnd using the embedded Runge-Kutta-Fehlberg 4(5) pair,
// accepting a step only when its local error estimate is within tol and
// otherwise shrinking h and retrying. Returns:
//   ts, ys — the accepted trajectory, one entry per accepted step plus t0
//   steps — array of { t0, t1, h, dy, deriv }, one entry per accepted step:
//     deriv is f(t0, y) at the step's start, dy is the accepted increment
//   acceptedCount, rejectedCount — step counters
rk45Solve = (f, t0, y0, tEnd, tol, maxSteps) => {
  const c2 = 1 / 4, c3 = 3 / 8, c4 = 12 / 13, c5 = 1, c6 = 1 / 2
  const a21 = 1 / 4
  const a31 = 3 / 32, a32 = 9 / 32
  const a41 = 1932 / 2197, a42 = -7200 / 2197, a43 = 7296 / 2197
  const a51 = 439 / 216, a52 = -8, a53 = 3680 / 513, a54 = -845 / 4104
  const a61 = -8 / 27, a62 = 2, a63 = -3544 / 2565, a64 = 1859 / 4104, a65 = -11 / 40
  const b1 = 25 / 216, b3 = 1408 / 2565, b4 = 2197 / 4104, b5 = -1 / 5
  const b1s = 16 / 135, b3s = 6656 / 12825, b4s = 28561 / 56430, b5s = -9 / 50, b6s = 2 / 55

  const direction = tEnd >= t0 ? 1 : -1
  const span = Math.abs(tEnd - t0)
  const hMin = span * 1e-10
  const hMax = span

  const ts = [t0]
  const ys = [y0]
  const steps = []
  let t = t0, y = y0
  let h = direction * span / 100
  let acceptedCount = 0
  let rejectedCount = 0

  // Safety valve: caps total solver attempts (accepted + rejected) so a
  // pathological (e.g. non-converging) right-hand side can't hang the
  // browser retrying rejected steps forever.
  const maxAttempts = maxSteps * 50
  let attempts = 0

  while (attempts < maxAttempts) {
    if (direction > 0 && t >= tEnd - 1e-12) break
    if (direction < 0 && t <= tEnd + 1e-12) break
    if (acceptedCount >= maxSteps) break
    attempts += 1

    if (direction > 0 && t + h > tEnd) h = tEnd - t
    if (direction < 0 && t + h < tEnd) h = tEnd - t

    const k1 = f(t, y)
    const k2 = f(t + c2 * h, y + h * (a21 * k1))
    const k3 = f(t + c3 * h, y + h * (a31 * k1 + a32 * k2))
    const k4 = f(t + c4 * h, y + h * (a41 * k1 + a42 * k2 + a43 * k3))
    const k5 = f(t + c5 * h, y + h * (a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4))
    const k6 = f(t + c6 * h, y + h * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5))

    const y4 = y + h * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5)
    const y5 = y + h * (b1s * k1 + b3s * k3 + b4s * k4 + b5s * k5 + b6s * k6)

    const error = Math.abs(y5 - y4)

    let factor = 5
    if (error !== 0) {
      factor = 0.9 * Math.pow(tol / error, 0.2)
      if (factor > 5) factor = 5
      if (factor < 0.1) factor = 0.1
    }
    let hNew = h * factor
    if (Math.abs(hNew) < hMin) hNew = direction * hMin
    if (Math.abs(hNew) > hMax) hNew = direction * hMax

    if (!Number.isFinite(error) || error <= tol || Math.abs(h) <= hMin) {
      const tNext = t + h
      steps.push({t0: t, t1: tNext, h, dy: y5 - y, deriv: k1})
      t = tNext
      y = y5
      ts.push(t)
      ys.push(y)
      acceptedCount += 1
      h = hNew
    } else {
      rejectedCount += 1
      h = hNew
    }
  }

  return {ts, ys, steps, acceptedCount, rejectedCount}
}

// result: runs adaptive RKF45 on [0, b] at the current tolerance ε.
//   steps, ts, ys — see rk45Solve
//   f — parsed right-hand side, reused for the reference trajectory in setup
//   acceptedCount, rejectedCount — step counters
//   finalValue — the computed y(b)
//   trueValue — reference y(b) via fixed-step RK4 with a very fine step count
//   actualError — |finalValue - trueValue|
result = {
  const math = window.math

  const f = VM.expressions.makeFunction2(math, committed.fText)
  const fallback = (t, y) => NaN

  if (!f) return {steps: [], ts: [], ys: [], f: fallback, acceptedCount: 0, rejectedCount: 0, finalValue: NaN, trueValue: NaN, actualError: NaN}

  const t0 = 0
  const y0 = committed.y0
  const b = committed.bEnd
  if (!Number.isFinite(y0) || !Number.isFinite(b) || b === t0) return {steps: [], ts: [], ys: [], f, acceptedCount: 0, rejectedCount: 0, finalValue: NaN, trueValue: NaN, actualError: NaN}

  const solved = rk45Solve(f, t0, y0, b, tolerance, maxSteps)

  if (solved.steps.length === 0) return {steps: [], ts: [], ys: [], f, acceptedCount: 0, rejectedCount: 0, finalValue: NaN, trueValue: NaN, actualError: NaN}

  const refSteps = 2000
  const trueValue = VM.numerical.rk4Solve(f, t0, y0, b, refSteps).ys[refSteps]
  const finalValue = solved.ys[solved.ys.length - 1]
  const actualError = Math.abs(finalValue - trueValue)

  return {steps: solved.steps, ts: solved.ts, ys: solved.ys, f, acceptedCount: solved.acceptedCount, rejectedCount: solved.rejectedCount, finalValue, trueValue, actualError}
}

// setup computes the data needed by the main chart: a fine reference
// trajectory for the smooth "true" solution curve, the adaptive trajectory,
// a per-point marker size (bigger where steps have shrunk, so refinement is
// visible directly on the solution curve), and the axis ranges.
setup = {
  if (result.steps.length === 0) {
    const empty = {lo: -3, hi: 3}
    return {ts: [], ys: [], pointSizes: [], refTs: [], refYs: [], tRange: empty, yRange: {lo: -1, hi: 1}}
  }

  const t0 = 0
  const y0 = committed.y0
  const b = committed.bEnd

  // A fine reference trajectory for the smooth "true" solution curve.
  const reference = VM.numerical.rk4Solve(result.f, t0, y0, b, 500)

  const tRange = VM.plotting.paddedRange([t0, b], {emptyRange: [-3, 3], relativePadding: 0.1, minPadding: 0.25})

  // Cap at 1e8 so a diverging solution doesn't blow up the axis range.
  const yValues = []
  for (const y of reference.ys) if (Number.isFinite(y) && Math.abs(y) <= 1e8) yValues.push(y)
  for (const y of result.ys) if (Number.isFinite(y) && Math.abs(y) <= 1e8) yValues.push(y)
  const yRange = VM.plotting.paddedRange(yValues, {emptyRange: [-1, 1], relativePadding: 0.25, minPadding: 0.5})

  // Marker size per accepted point, on a log scale of the step size that
  // produced it — smaller steps (more refinement) get bigger markers. t0
  // has no incoming step, so it reuses the first step's size.
  let minH = Infinity
  let maxH = 0
  for (const step of result.steps) {
    const h = Math.abs(step.h)
    if (h < minH) minH = h
    if (h > maxH) maxH = h
  }
  if (minH === maxH) {
    minH = minH * 0.5
    maxH = maxH * 2
  }
  const logMin = Math.log10(minH)
  const logMax = Math.log10(maxH)

  const sizeForH = h => {
    let ratio = (logMax - Math.log10(h)) / (logMax - logMin)
    if (!Number.isFinite(ratio)) ratio = 0.5
    return 5 + 9 * ratio
  }

  const pointSizes = [sizeForH(Math.abs(result.steps[0].h))]
  for (const step of result.steps) pointSizes.push(sizeForH(Math.abs(step.h)))

  return {ts: result.ts, ys: result.ys, pointSizes, refTs: reference.ts, refYs: reference.ys, tRange, yRange}
}

// mainPlot draws the fine reference solution curve and the adaptive RK45
// trajectory, with marker size showing where the step size shrunk to
// resolve fast dynamics.
//
// _plotDiv is kept in a closure so Plotly reuses the same DOM node across
// reactive updates — this preserves zoom and pan when the tolerance or max
// steps changes.
mainPlot = {
  const Plotly = window.Plotly
  let _plotDiv = null

  return (data, showRefSolution, showAdaptive) => {
    const traces = []

    if (showRefSolution) {
      traces.push({x: data.refTs, y: data.refYs, type: "scatter", mode: "lines", name: "Reference solution", showlegend: false, line: {color: chartColors.fn, width: 3}})
    }

    if (showAdaptive && data.ts.length > 0) {
      traces.push({
        x: data.ts, y: data.ys, type: "scatter", mode: "lines+markers", name: "Adaptive RK45", showlegend: false,
        line: {color: chartColors.alt, width: 2}, marker: {color: chartColors.alt, size: data.pointSizes, line: {color: chartColors.halo, width: 1}}
      })
    }

    const layout = {
      xaxis: {title: "t", range: [data.tRange.lo, data.tRange.hi], zeroline: true},
      yaxis: {title: "y", range: [data.yRange.lo, data.yRange.hi], zeroline: true},
      hovermode: "closest",
      uirevision: "static",  // keeps zoom/pan state across reactive updates
      autosize: true
    }
    const config = VM.plotting.config()

    // First render: create the div. Later renders: update in place (preserves zoom).
    if (!_plotDiv) {
      _plotDiv = document.createElement("div")
      _plotDiv.className = "plotly-box-large"
      Plotly.newPlot(_plotDiv, traces, layout, config)
      // newPlot fires before the div is in the DOM, so Plotly measures 0 width.
      // A ResizeObserver catches the real size once it's inserted and laid out,
      // and keeps correcting it on any later layout/viewport change.
      VM.plotting.autoResize(_plotDiv)
    } else {
      Plotly.react(_plotDiv, traces, layout, config)
    }
    return _plotDiv
  }
}

// maxSteps parses the Max steps field's expression, falling back to 2000
// for invalid input so the solver always gets a usable step limit.
maxSteps = {
  const math = window.math
  const parsed = VM.expressions.makeNumber(math, maxStepsInput)
  if (parsed === null || !Number.isFinite(parsed)) return 2000
  return Math.max(1, Math.round(parsed))
}

mainPlot(setup, showRefSolution, showAdaptive)

viewof methods = {
  const items = []
  for (const label of methodOptions) items.push({label, color: methodColors.get(label)})
  return VM.ui.legendOverlay(items, {value: methodOptions})
}

viewof tolExponent = {
  const tolSlider = Inputs.range([0, 10], {step: 0.5, value: 0, label: "Tolerance exponent e (ε = 10⁻ᵉ)"})
  tolSlider.classList.add("ojs-fill")
  return tolSlider
}

Step table

// Step table — one row per accepted step, left to right.
iterationTable = {
  if (result.steps.length === 0) {
    // Always return a Node so Quarto's OJS runtime never latches this
    // declaration cell as hidden.
    return document.createElement("div")
  }

  const formattedRows = []
  for (let i = 0; i < result.steps.length; i++) {
    const step = result.steps[i]
    formattedRows.push([
      i,
      step.t0.toPrecision(10),
      step.t1.toPrecision(10),
      step.h.toExponential(4),
      step.dy.toExponential(4),
      step.deriv.toExponential(4)
    ])
  }

  return VM.ui.renderTable({
    html,
    headers: ["i", tex`t_i`, tex`t_{i+1}`, tex`h_i`, tex`\Delta y_i`, tex`y'(t_i)`],
    csvHeaders: ["i", "t_i", "t_{i+1}", "h_i", "dy_i", "yprime_t_i"],
    filename: "adaptive-rk45-steps.csv",
    rows: formattedRows
  })
}

How it works

The explicit methods take every step with the same size \(h\). An adaptive method instead picks a new \(h\) at each step: small where the solution is changing quickly, large where it is nearly straight.

To decide, it needs an estimate of the error a step makes. RKF45 (Runge–Kutta–Fehlberg) gets one cheaply: from the same six slope evaluations it forms two approximations of \(y(t + h)\), one of 4th order and one of 5th order. Their difference

\[ E = \bigl|\, y_{\text{5th}} - y_{\text{4th}} \,\bigr| \]

is the error estimate for that step. It is then compared with the tolerance \(\varepsilon\) set by the slider:

Accept if \(E \le \varepsilon\): advance to \(t + h\) using the 5th-order value.
Reject if \(E > \varepsilon\): stay at \(t\) and retry with a smaller \(h\).

Either way the next step size is \(h \cdot 0.9\,(\varepsilon / E)^{1/5}\), so a step well inside tolerance grows \(h\) and a step outside it shrinks \(h\). The factor is capped between \(0.1\) and \(5\) so a single step can’t change \(h\) too drastically.

A larger \(\varepsilon\) means fewer, bigger steps and a rougher curve; a smaller \(\varepsilon\) means more steps and a curve that hugs the reference solution.