Adaptive Integration

Refines a numerical integral by subdividing only the sub-intervals where the estimate is still inaccurate.

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.numerical.lagrangeQuadratic, VM.expressions.makeFunction, VM.expressions.makeNumber, VM.plotting.paddedRange, VM.ui.renderTable, VM.numerical.simpsonEstimate, VM.ui.applyExampleParams.
VM = window.VM

traceOptions = ["f(x)", "Adaptive Simpson panels", "Partition points"]

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

traceColors = new Map([
  ["f(x)", chartColors.fn],
  ["Adaptive Simpson panels", chartColors.ok],
  ["Partition points", chartColors.ink]
])

// examples backs the "Try an example" dropdown near the top controls.
examples = [
  {
    title: "Oscillation near an endpoint",
    params: {max: "15", f: "sin(1 / (x + 0.3))", a: "0", b: "pi"}
  },
  {
    title: "Challenging oscillations",
    params: {max: "15", f: "sin(1/x)", a: "0.05", b: "1"}
  },
  {
    title: "Narrow spike",
    params: {max: "15", f: "1/(0.01 + x^2)", a: "-1", b: "1"}
  },
  {
    title: "Steep transition",
    params: {max: "15", f: "tanh(20*(x - 0.2))", a: "-1", b: "1"}
  },
  {
    title: "Nonsmooth corner",
    params: {max: "15", f: "abs(x)", a: "-2", b: "2"}
  },
  {
    title: "Endpoint behavior",
    params: {max: "15", f: "sqrt(x)", a: "0", b: "1"}
  }
]

viewof fText = {
  const params = new URLSearchParams(window.location.search)
  const fTextInput = Inputs.text({label: "Function f(x)", value: params.get("f") ?? "sin(1 / (x + 0.3))", placeholder: "Example: sin(x), x^3 - x, 1 / (1 + x^2)"})
  fTextInput.classList.add("ojs-fill")
  fTextInput.dataset.exampleField = "f"
  return fTextInput
}

// Plot commits the current function/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 a0 = {
  const params = new URLSearchParams(window.location.search)
  const a0Input = Inputs.text({label: "Left endpoint a", value: params.get("a") ?? "0", placeholder: "Example: 0, -pi"})
  a0Input.dataset.exampleField = "a"
  return a0Input
}

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

// Always-live (not gated behind Plot) — it bounds the recursion depth
// 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 maxDepthInput = {
  const input = Inputs.text({value: initialMax, label: "Max depth", placeholder: "Example: 15, 5*3"})
  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"]',
  a: '[data-example-field="a"]',
  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/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 a0, viewof b0]) {
    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/a0/b0 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 a0`/`viewof b0` (the stable DOM views) rather than
// `fText`/`a0`/`b0` (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 aText = (viewof a0).value
  const bText = (viewof b0).value
  const aVal = VM.expressions.makeNumber(math, aText)
  const bVal = VM.expressions.makeNumber(math, bText)

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

  return {fText: fVal, a0: aVal, b0: 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", maxDepthInput)
  history.replaceState(null, "", `${window.location.pathname}?${params}`)
  return true
}

// result: runs adaptive Simpson quadrature on [a, b] at the current
// tolerance ε (recursing until each leaf's local error estimate is within
// its share of ε, or maxDepth is reached).
//   leaves — array of { a, b, mid, fa, fmid, fb, estimate, errorEstimate,
//     depth, width }, one entry per accepted subinterval, left to right
//   f — parsed function, reused for curve sampling in setup
//   totalEstimate — sum of leaf estimates (the adaptive quadrature's answer)
//   totalErrorEstimate — sum of leaf Richardson error estimates
//   trueIntegral — reference integral via a very fine composite Simpson's rule
//   actualError — |totalEstimate - trueIntegral|
result = {
  const math = window.math

  const f = VM.expressions.makeFunction(math, committed.fText)
  const fallback = x => NaN

  if (!f) return {leaves: [], f: fallback, totalEstimate: NaN, totalErrorEstimate: NaN, trueIntegral: NaN, actualError: NaN}

  const a = committed.a0
  const b = committed.b0
  if (!Number.isFinite(a) || !Number.isFinite(b) || a === b) return {leaves: [], f, totalEstimate: NaN, totalErrorEstimate: NaN, trueIntegral: NaN, actualError: NaN}

  const lo = Math.min(a, b)
  const hi = Math.max(a, b)

  const fa = f(lo)
  const fb = f(hi)
  if (!Number.isFinite(fa) || !Number.isFinite(fb)) return {leaves: [], f, totalEstimate: NaN, totalErrorEstimate: NaN, trueIntegral: NaN, actualError: NaN}

  const depthLimit = maxDepth
  // Safety valve: caps total leaves so a pathological (e.g. non-converging
  // or discontinuous) function can't hang the browser with an exponential
  // blow-up of subdivisions.
  const maxLeaves = 4000

  // Simpson's rule on [x0, x1], reusing the already-evaluated endpoint
  // values y0, y1 and computing only the new midpoint value.
  const simpsonPiece = (x0, x1, y0, y1) => {
    const mid = (x0 + x1) / 2
    const ymid = f(mid)
    const value = (x1 - x0) / 6 * (y0 + 4 * ymid + y1)
    return {mid, ymid, value}
  }

  const leaves = []

  const recurse = (a0, b0, fa0, fmid0, fb0, whole, tol, depth) => {
    const c = (a0 + b0) / 2

    if (leaves.length >= maxLeaves) {
      leaves.push({a: a0, b: b0, mid: c, fa: fa0, fmid: fmid0, fb: fb0, estimate: whole, errorEstimate: 0, depth, width: b0 - a0})
      return
    }

    const left = simpsonPiece(a0, c, fa0, fmid0)
    const right = simpsonPiece(c, b0, fmid0, fb0)
    const combined = left.value + right.value
    const errorEst = (combined - whole) / 15

    if (!Number.isFinite(combined) || depth >= depthLimit || Math.abs(combined - whole) <= 15 * tol) {
      leaves.push({a: a0, b: b0, mid: c, fa: fa0, fmid: fmid0, fb: fb0, estimate: combined + errorEst, errorEstimate: Math.abs(errorEst), depth, width: b0 - a0})
      return
    }

    recurse(a0, c, fa0, left.ymid, fmid0, left.value, tol / 2, depth + 1)
    recurse(c, b0, fmid0, right.ymid, fb0, right.value, tol / 2, depth + 1)
  }

  const whole0 = simpsonPiece(lo, hi, fa, fb)
  recurse(lo, hi, fa, whole0.ymid, fb, whole0.value, tolerance, 0)

  let totalEstimate = 0
  let totalErrorEstimate = 0
  for (const leaf of leaves) {
    totalEstimate += leaf.estimate
    totalErrorEstimate += leaf.errorEstimate
  }

  const trueIntegral = VM.numerical.simpsonEstimate(f, lo, hi, 4000)
  const actualError = Math.abs(totalEstimate - trueIntegral)

  return {leaves, f, totalEstimate, totalErrorEstimate, trueIntegral, actualError}
}

// setup computes the data needed by the main chart: a dense quadratic
// (Simpson) sample per leaf for shading, the partition boundary points, and
// the axis ranges / smooth f(x) curve.
setup = {
  if (result.leaves.length === 0) {
    const empty = {lo: -3, hi: 3}
    return {panels: [], nodesX: [], nodesY: [], boundaries: [], maxDepthUsed: 0, sampleXs: [], sampleYs: [], xRange: empty, xBuffer: empty, yRange: {lo: -1, hi: 1}}
  }

  let maxDepthUsed = 0
  for (const leaf of result.leaves) {
    if (leaf.depth > maxDepthUsed) maxDepthUsed = leaf.depth
  }

  // Dense parabola samples per leaf, via the Lagrange quadratic through its
  // (a, fa), (mid, fmid), (b, fb) — the same curve Simpson's rule integrates
  // exactly on that leaf.
  const samplesPerLeaf = 16
  const panels = []
  for (const leaf of result.leaves) {
    const sample = VM.numerical.lagrangeQuadratic(leaf.a, leaf.fa, leaf.mid, leaf.fmid, leaf.b, leaf.fb, samplesPerLeaf)
    panels.push({xs: sample.xs, ys: sample.ys, depth: leaf.depth})
  }

  // Partition boundary points — every leaf's left edge, plus the final
  // right edge — used for the node markers and the vertical divider lines.
  const nodesX = [result.leaves[0].a]
  const nodesY = [result.leaves[0].fa]
  for (const leaf of result.leaves) {
    nodesX.push(leaf.b)
    nodesY.push(leaf.fb)
  }
  const boundaries = nodesX.slice(1, -1)

  const lo = result.leaves[0].a
  const hi = result.leaves[result.leaves.length - 1].b
  const xRange = VM.plotting.paddedRange([lo, hi], {emptyRange: [-3, 3], relativePadding: 0.1, minPadding: 0.25})

  // A wider domain (double the padding of xRange) so the curve is already
  // sampled beyond the initial view — panning/zooming out reveals more of
  // the curve instead of blank canvas. The initial view stays xRange.
  const xBuffer = VM.plotting.paddedRange([lo, hi], {emptyRange: [-6, 6], relativePadding: 0.25, minPadding: 0.5})

  // Sample f across the buffer range for the smooth curve trace.
  const sampleCount = 700
  const sampleXs = Array.from({length: sampleCount}, (_, i) => xBuffer.lo + (xBuffer.hi - xBuffer.lo) * i / (sampleCount - 1))
  const sampleYs = []
  for (const x of sampleXs) {
    const y = result.f(x)
    if (Number.isFinite(y) && Math.abs(y) < 1e8) {
      sampleYs.push(y)
    } else {
      sampleYs.push(null)
    }
  }

  // y-range from f over [lo, hi] plus node and panel values.
  const yValues = [0]
  for (let i = 0; i < sampleXs.length; i++) {
    const x = sampleXs[i]
    if (x >= xRange.lo && x <= xRange.hi && sampleYs[i] !== null) yValues.push(sampleYs[i])
  }
  for (const y of nodesY) if (Number.isFinite(y)) yValues.push(y)
  for (const panel of panels) {
    for (const y of panel.ys) if (Number.isFinite(y)) yValues.push(y)
  }
  const yRange = VM.plotting.paddedRange(yValues, {emptyRange: [-1, 1], relativePadding: 0.25, minPadding: 0.5})

  return {panels, nodesX, nodesY, boundaries, maxDepthUsed, sampleXs, sampleYs, xRange, xBuffer, yRange}
}

// mainPlot draws f(x) and, for each leaf of the adaptive partition, the
// quadratic panel Simpson's rule integrates exactly there — shaded more
// darkly the deeper that leaf was subdivided, so heavily refined regions
// (where f is hardest to approximate) stand out. Vertical dividers and dot
// markers show the partition points themselves.
//
// _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
// depth changes.
mainPlot = {
  const Plotly = window.Plotly
  let _plotDiv = null

  return (data, showFCurve, showPanels, showPartitionPoints) => {
    const traces = []

    if (showFCurve) {
      // The function curve
      traces.push({x: data.sampleXs, y: data.sampleYs, type: "scatter", mode: "lines", name: "f(x)", showlegend: false, line: {color: chartColors.fn, width: 3}})
    }

    if (showPanels) {
      for (let i = 0; i < data.panels.length; i++) {
        const panel = data.panels[i]
        let depthRatio = 0
        if (data.maxDepthUsed > 0) depthRatio = panel.depth / data.maxDepthUsed
        const fillAlpha = 0.10 + 0.45 * depthRatio

        traces.push({
          x: panel.xs, y: panel.ys, type: "scatter", mode: "lines", name: "Adaptive Simpson panels", showlegend: false,
          fill: "tozeroy", fillcolor: VM.plotting.alpha("ok", fillAlpha), line: {color: chartColors.ok, width: 1}
        })
      }
    }

    if (showPartitionPoints && data.nodesX.length > 0) {
      traces.push({x: data.nodesX, y: data.nodesY, type: "scatter", mode: "markers", name: "Partition points", showlegend: false, marker: {color: chartColors.ink, size: 6, symbol: "circle", line: {color: chartColors.halo, width: 1}}})
    }

    // Vertical dividers at interior partition points, spanning the full
    // plot height, so unevenly sized panels are easy to see. Tied to the
    // Partition points toggle, matching how newton-method tied its dotted
    // guide lines to showNewtonPoints.
    const shapes = []
    if (showPartitionPoints) {
      for (const x of data.boundaries) {
        shapes.push({type: "line", xref: "x", yref: "paper", x0: x, x1: x, y0: 0, y1: 1, line: {color: "rgba(17, 24, 39, 0.25)", width: 1, dash: "dot"}})
      }
    }

    const layout = {
      xaxis: {title: "x", range: [data.xRange.lo, data.xRange.hi], zeroline: true},
      yaxis: {title: "y", range: [data.yRange.lo, data.yRange.hi], zeroline: true},
      shapes,
      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
  }
}

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

mainPlot(setup, showFCurve, showPanels, showPartitionPoints)

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

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
}

Partition table

// Partition table — one row per leaf of the accepted adaptive partition,
// left to right.
iterationTable = {
  if (result.leaves.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.leaves.length; i++) {
    const leaf = result.leaves[i]
    formattedRows.push([
      i,
      leaf.a.toPrecision(10),
      leaf.b.toPrecision(10),
      leaf.width.toExponential(4),
      leaf.depth,
      leaf.estimate.toPrecision(10),
      leaf.errorEstimate.toExponential(4)
    ])
  }

  return VM.ui.renderTable({
    html,
    headers: ["i", tex`a_i`, tex`b_i`, "width", "depth", tex`S_i`, tex`|E_i|`],
    csvHeaders: ["i", "a_i", "b_i", "width", "depth", "S_i (leaf estimate)", "|E_i| (leaf error estimate)"],
    filename: "adaptive-integration-partition.csv",
    rows: formattedRows
  })
}

How it works

Numerical integration lays a fixed grid over \([a,b]\) and uses the same panel width everywhere. An adaptive method instead subdivides only where it has to: narrow panels where \(f\) wiggles or turns sharply, wide ones where \(f\) is nearly a parabola already.

To decide, it needs an estimate of the error a panel makes. Simpson’s rule on \([p,q]\) with midpoint \(m = (p+q)/2\) is

\[ S[p,q] = \frac{q-p}{6}\Bigl(f(p) + 4f(m) + f(q)\Bigr), \]

and the same panel split in half gives a second, finer estimate \(S[p,m] + S[m,q]\). Since Simpson’s error scales like \(h^4\), halving \(h\) shrinks it by a factor of \(16\), so the difference between the two is almost all error:

\[ E = \frac{\bigl|\, S[p,m] + S[m,q] - S[p,q] \,\bigr|}{15}. \]

It is then compared with the tolerance \(\varepsilon\) set by the slider:

Accept if \(E \le \varepsilon\): keep \([p,q]\) as a panel of the final partition, using the finer estimate.
Refine if \(E > \varepsilon\): bisect at \(m\) and recurse on each half with \(\varepsilon/2\), so the total error stays bounded by the original tolerance.

A larger \(\varepsilon\) means fewer, wider panels and a coarser partition; a smaller \(\varepsilon\) means more panels, clustered exactly where \(f\) is hardest to approximate.