Adaptive Methods for ODEs

// 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 = window.VM

The classical RK4 method in Explicit Methods advances with a fixed step size \(h = b/n\) everywhere — spending the same effort even where the solution is nearly linear. RKF45 (Runge–Kutta–Fehlberg) instead adapts \(h\) at every step, driven by a target tolerance \(\varepsilon\) rather than a fixed step count.

From the same six slope evaluations,

\[ \begin{aligned} k_1 &= f(t_i, y_i), \\ k_2 &= f\!\left(t_i + \tfrac{1}{4}h,\ y_i + \tfrac{1}{4}h k_1\right), \\ k_3 &= f\!\left(t_i + \tfrac{3}{8}h,\ y_i + h\left(\tfrac{3}{32}k_1 + \tfrac{9}{32}k_2\right)\right), \\ k_4 &= f\!\left(t_i + \tfrac{12}{13}h,\ y_i + h\left(\tfrac{1932}{2197}k_1 - \tfrac{7200}{2197}k_2 + \tfrac{7296}{2197}k_3\right)\right), \\ k_5 &= f\!\left(t_i + h,\ y_i + h\left(\tfrac{439}{216}k_1 - 8k_2 + \tfrac{3680}{513}k_3 - \tfrac{845}{4104}k_4\right)\right), \\ k_6 &= f\!\left(t_i + \tfrac{1}{2}h,\ y_i + h\left(-\tfrac{8}{27}k_1 + 2k_2 - \tfrac{3544}{2565}k_3 + \tfrac{1859}{4104}k_4 - \tfrac{11}{40}k_5\right)\right), \end{aligned} \]

RKF45 builds two embedded estimates of \(y_{i+1}\) — a fourth-order one and a fifth-order one:

\[ y_{i+1}^{(4)} = y_i + h\left(\tfrac{25}{216}k_1 + \tfrac{1408}{2565}k_3 + \tfrac{2197}{4104}k_4 - \tfrac{1}{5}k_5\right), \]

\[ y_{i+1}^{(5)} = y_i + h\left(\tfrac{16}{135}k_1 + \tfrac{6656}{12825}k_3 + \tfrac{28561}{56430}k_4 - \tfrac{9}{50}k_5 + \tfrac{2}{55}k_6\right). \]

Their difference \(E = \left|y_{i+1}^{(5)} - y_{i+1}^{(4)}\right|\) estimates the local error. If \(E \le \varepsilon\), the step is accepted using \(y_{i+1}^{(5)}\) (the higher-order “local extrapolation”), and the next step size grows by a factor of roughly \((\varepsilon/E)^{1/5}\). If \(E > \varepsilon\), the step is rejected and retried with a shrunk \(h\) — so every accepted step satisfies the same local error tolerance \(\varepsilon\), no matter how the dynamics vary across \([0, b]\).

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

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

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

// 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 and max steps
// field below are unaffected and stay reactive.
viewof plotTrigger = {
  const plotButton = Inputs.button("Plot", {value: 0, reduce: v => v + 1})
  plotButton.classList.add("ojs-auto")
  return plotButton
}

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

viewof maxStepsInput = {
  const input = Inputs.text({value: "2000", label: "Max steps", placeholder: "Example: 2000, 1000*2"})
  input.classList.add("ojs-auto")
  return input
}

// 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))
}

tolerance = Math.pow(10, -tolExponent)

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

// 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) => {
    const traces = [
      // The fine reference solution curve
      {x: data.refTs, y: data.refYs, type: "scatter", mode: "lines", name: "Reference solution", showlegend: true, line: {color: "#2563eb", width: 3}}
    ]

    if (data.ts.length > 0) {
      traces.push({
        x: data.ts, y: data.ys, type: "scatter", mode: "lines+markers", name: "Adaptive RK45", showlegend: true,
        line: {color: "#dc2626", width: 2}, marker: {color: "#dc2626", size: data.pointSizes, line: {color: "white", width: 1}}
      })
    }

    const layout = {
      margin: {l: 0, r: 0, t: 0, b: 0},
      xaxis: {title: "t", range: [data.tRange.lo, data.tRange.hi], zeroline: true},
      yaxis: {title: "y", range: [data.yRange.lo, data.yRange.hi], zeroline: true},
      legend: {orientation: "h", y: -0.08},
      hovermode: "closest",
      uirevision: "static",  // keeps zoom/pan state across reactive updates
      autosize: true
    }
    const config = {responsive: true, displaylogo: false}

    // 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.
      new ResizeObserver(() => Plotly.Plots.resize(_plotDiv)).observe(_plotDiv)
    } else {
      Plotly.react(_plotDiv, traces, layout, config)
    }
    return _plotDiv
  }
}

// summary reports the headline numbers for the current tolerance: how many
// steps were accepted and rejected, the resulting y(b), and the actual
// error against a fine reference solution.
//
// Defined (and displayed) right here, rather than as a hidden cell invoked
// from elsewhere: it already returns the final displayable node, and a
// `//| output: false` cell still attaches its returned node to its own
// (hidden) output slot — a later bare reference to it would then find the
// node already parented elsewhere and fall back to showing an inspector
// summary instead of the node itself.
summary = {
  const container = document.createElement("p")

  if (result.steps.length === 0) {
    container.textContent = "Enter a valid y', initial condition, and endpoint, then click Plot."
    return container
  }

  const stats = [
    ["Tolerance ε", tolerance.toExponential(2)],
    ["Steps accepted", String(result.acceptedCount)],
    ["Steps rejected", String(result.rejectedCount)],
    ["y(b) estimate", result.finalValue.toPrecision(8)],
    ["Actual error", result.actualError.toExponential(3)]
  ]

  for (let i = 0; i < stats.length; i++) {
    const [label, value] = stats[i]
    if (i > 0) container.append(document.createTextNode("  ·  "))
    const strong = document.createElement("strong")
    strong.textContent = `${label}: `
    container.append(strong, document.createTextNode(value))
  }

  return container
}

mainPlot(setup)

The chart below shows each accepted step’s increment \(\Delta y_i = y_{i+1}^{(5)} - y_i\) across its time interval \([t_i, t_{i+1}]\), with a semi-transparent overlay of the derivative \(y'(t) = f(t, y(t))\) evaluated along the solution. The two live on different scales — \(\Delta y_i \approx h_i\, y'(t_i)\), and \(h_i\) is usually well below 1 — so the derivative curve is rescaled to the increments’ own amplitude before overlaying; only its shape, not its absolute value, should be compared against \(\Delta y_i\) (hover it for the true value).

// incrementsPlot draws the accepted increment Δy_i as a horizontal segment
// across each step's time interval, with a semi-transparent filled overlay
// of the derivative y'(t) = f(t, y(t)) evaluated along the adaptive
// trajectory's own points (so it shares the same time resolution as the
// increments next to it — sparser where steps are large, denser where
// they've shrunk).
//
// Δy and y' generally live on very different scales (Δy ≈ h·y', and h is
// usually « 1), so plotting them raw on shared axes lets whichever one is
// larger swamp the other. The derivative curve is rescaled to match the
// increments' own peak amplitude before overlaying — only its shape (where
// it rises/falls/crosses zero), not its absolute height, is meant to be
// compared against Δy. The true (unscaled) derivative value is kept for the
// tooltip.
//
// Defined (and displayed) right here for the same reason as `summary`
// above — it already returns the final displayable node.
incrementsPlot = {
  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")
  }

  let maxAbsDy = 0
  for (const step of result.steps) {
    const absDy = Math.abs(step.dy)
    if (absDy > maxAbsDy) maxAbsDy = absDy
  }

  let maxAbsDeriv = 0
  const rawDeriv = []
  for (let i = 0; i < result.ts.length; i++) {
    const value = result.f(result.ts[i], result.ys[i])
    rawDeriv.push(value)
    if (Number.isFinite(value) && Math.abs(value) > maxAbsDeriv) maxAbsDeriv = Math.abs(value)
  }

  // Scale the derivative so its peak matches the increments' peak — a pure
  // shape-matching overlay, not a shared-unit axis.
  let derivScale = 1
  if (maxAbsDeriv > 0) derivScale = maxAbsDy / maxAbsDeriv

  const derivPoints = []
  for (let i = 0; i < result.ts.length; i++) {
    derivPoints.push({t: result.ts[i], trueY: rawDeriv[i], y: rawDeriv[i] * derivScale})
  }

  const incrementName = "Increment Δy"
  const derivName = "y′(t) (derivative, rescaled to match)"

  const tooltipOf = step => `[${step.t0.toPrecision(6)}, ${step.t1.toPrecision(6)}]  h=${step.h.toExponential(3)}  Δy=${step.dy.toExponential(3)}`
  const derivTooltipOf = d => `t=${d.t.toPrecision(6)}  y′(t)=${d.trueY.toExponential(3)}`

  const colorDomain = [derivName, incrementName]
  const colorRange = ["#2563eb", "#dc2626"]

  const chart = Plot.plot({
    width,
    height: 320,
    marginLeft: 46,
    marginBottom: 30,
    x: {label: "t"},
    y: {label: "Δy  (y′ rescaled to match)"},
    color: {domain: colorDomain, range: colorRange},
    marks: [
      Plot.areaY(derivPoints, {x: "t", y: "y", fill: () => derivName, fillOpacity: 0.12, title: derivTooltipOf}),
      Plot.lineY(derivPoints, {x: "t", y: "y", stroke: () => derivName, strokeOpacity: 0.3, strokeWidth: 1.5, title: derivTooltipOf}),
      Plot.ruleY([0], {stroke: "#94a3b8", strokeWidth: 1}),
      Plot.ruleY(result.steps, {
        x1: "t0", x2: "t1", y: "dy",
        stroke: () => incrementName,
        strokeWidth: 5,
        strokeLinecap: "round",
        title: tooltipOf
      })
    ]
  })
  const legend = Plot.legend({color: {domain: colorDomain, range: colorRange}})

  const container = document.createElement("div")
  container.append(chart, legend)
  return container
}

NoteStep 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
  })
}