Numerical Integration

Compares Riemann sums, the trapezoid rule, and Simpson’s rule for approximating a definite integral.

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

urlParams = new URLSearchParams(window.location.search)

methodOptions = ["f(x)", "Left Riemann sum", "Midpoint Riemann sum", "Trapezoidal rule", "Simpson's rule", "Interpolation nodes"]

methodSlugToLabel = new Map([
  ["f", "f(x)"],
  ["fx", "f(x)"],
  ["curve", "f(x)"],
  ["function", "f(x)"],

  ["left", "Left Riemann sum"],
  ["left-riemann", "Left Riemann sum"],
  ["left-riemann-sum", "Left Riemann sum"],
  ["left_riemann", "Left Riemann sum"],
  ["left_riemann_sum", "Left Riemann sum"],

  ["midpoint", "Midpoint Riemann sum"],
  ["mid", "Midpoint Riemann sum"],
  ["midpoint-riemann", "Midpoint Riemann sum"],
  ["midpoint-riemann-sum", "Midpoint Riemann sum"],
  ["midpoint_riemann", "Midpoint Riemann sum"],
  ["midpoint_riemann_sum", "Midpoint Riemann sum"],

  ["trapezoidal", "Trapezoidal rule"],
  ["trapezoid", "Trapezoidal rule"],
  ["trap", "Trapezoidal rule"],
  ["trapezoidal-rule", "Trapezoidal rule"],
  ["trapezoidal_rule", "Trapezoidal rule"],

  ["simpson", "Simpson's rule"],
  ["simpsons", "Simpson's rule"],
  ["simpson-rule", "Simpson's rule"],
  ["simpsons-rule", "Simpson's rule"],
  ["simpson_rule", "Simpson's rule"],
  ["simpsons_rule", "Simpson's rule"],

  ["nodes", "Interpolation nodes"],
  ["interpolation-nodes", "Interpolation nodes"],
  ["interpolation_nodes", "Interpolation nodes"]
])

methodLabelToSlug = new Map([
  ["f(x)", "f"],
  ["Left Riemann sum", "left"],
  ["Midpoint Riemann sum", "midpoint"],
  ["Trapezoidal rule", "trapezoidal"],
  ["Simpson's rule", "simpson"],
  ["Interpolation nodes", "nodes"]
])

// Matches the fill/line colors mainPlot uses for each method's trace, so
// the overlay legend's swatches double as a color key. #8b5cf6 (violet) has
// no VM.plotting.colors() token -- left as a literal.
// 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([
  ["f(x)", chartColors.fn],
  ["Left Riemann sum", "#8b5cf6"],
  ["Midpoint Riemann sum", chartColors.warn],
  ["Trapezoidal rule", chartColors.alt],
  ["Simpson's rule", chartColors.ok],
  ["Interpolation nodes", chartColors.ink]
])

initialMethods = {
  const raw = urlParams.get("show") ?? urlParams.get("methods")
  if (raw === null || raw.trim() === "") return methodOptions

  const selected = raw
    .split(",")
    .map(s => s.trim().toLowerCase())
    .map(s => methodSlugToLabel.get(s))
    .filter(Boolean)

  return selected.length > 0 ? Array.from(new Set(selected)) : methodOptions
}

initialN = {
  const raw = Number.parseInt(urlParams.get("n"), 10)
  return Number.isFinite(raw) ? Math.max(1, raw) : 6
}

initialMax = urlParams.get("max") ?? "40"

// examples backs the "Try an example" dropdown near the top controls.
//
// Each params object lists "max" before "n" — nControl's range depends
// reactively on maxN, so its view gets recreated when max changes;
// VM.ui.applyExampleParams applies params in this order with a yield
// between each, so the freshly-recreated slider is what "n" ends up
// applied to.
examples = [
  {
    title: "Signed area cancellation",
    params: {f: "x^3 - x", a: "-2", b: "2", max: "40", n: "6"}
  },
  {
    title: "Oscillatory integral",
    params: {f: "sin(5*x)", a: "0", b: "2*pi", max: "80", n: "12"}
  },
  {
    title: "Nonsmooth corner",
    params: {f: "abs(x)", a: "-2", b: "2", max: "50", n: "8"}
  },
  {
    title: "Rational peak",
    params: {f: "1/(1 + x^2)", a: "-5", b: "5", max: "60", n: "10"}
  },
  {
    title: "Sine area on one half-period",
    params: {f: "sin(x)", a: "0", b: "pi", max: "40", n: "6"}
  },
  {
    title: "Simpson's rule needs even n",
    params: {f: "sin(x)", a: "0", b: "pi", max: "35", n: "5"}
  }
]

viewof fText = {
  const params = new URLSearchParams(window.location.search)
  const fTextInput = Inputs.text({label: "Function f(x)", value: params.get("f") ?? "x^3 - x", 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 n slider and the method legend overlaying the chart
// 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 a0 = {
  const params = new URLSearchParams(window.location.search)
  const a0Input = Inputs.text({label: "Left endpoint a", value: params.get("a") ?? "-2", 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") ?? "2", placeholder: "Example: pi, 3.14"})
  b0Input.dataset.exampleField = "b"
  return b0Input
}

// Always-live (not gated behind Plot) — it bounds the step slider 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 maxNInput = {
  const input = Inputs.text({value: initialMax, label: "Max n", placeholder: "Example: 20, 5*4"})
  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 (nControl's view is reactively
// recreated whenever maxN changes, so referencing it by name here would
// make this cell -- and therefore applyExampleFromSelect -- reactive on
// it, forming a feedback loop with the fact that applying an example is
// what changes maxN in the first place).
exampleFieldSelectors = ({
  f: '[data-example-field="f"]',
  a: '[data-example-field="a"]',
  b: '[data-example-field="b"]',
  max: '[data-example-field="max"]',
  n: '[data-example-field="n"]'
})

// Applies the selected example's params to the fields above (including the
// n/max fields overlaid on the chart below) 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)
  })
}

showFCurve = methods.includes("f(x)")

showNodes = methods.includes("Interpolation nodes")

showLeftRiemann = methods.includes("Left Riemann sum")

showMidpoint = methods.includes("Midpoint Riemann sum")

showTrapezoidal = methods.includes("Trapezoidal rule")

showSimpson = methods.includes("Simpson's rule")

// Keeps n, max, and method choices shareable as the reactive controls change.
urlSyncControls = {
  const params = new URLSearchParams(window.location.search)
  params.set("n", String(nControl))
  params.set("max", maxNInput)
  params.set("show", methods.map(m => methodLabelToSlug.get(m)).filter(Boolean).join(","))
  history.replaceState(null, "", `${window.location.pathname}?${params}`)
  return true
}

// 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)
  params.set("n", String(nControl))
  params.set("max", maxNInput)
  params.set("show", methods.map(m => methodLabelToSlug.get(m)).filter(Boolean).join(","))
  history.replaceState(null, "", `${window.location.pathname}?${params}`)

  return {fText: fVal, a0: aVal, b0: bVal}
}

// result: for the current committed function/endpoints, builds the
// quadrature history used by the convergence plots and table.
//   rows — array of { n, leftArea, midpointArea, trapezoidalArea, simpsonArea,
//     trueIntegral, leftError, midpointError, trapezoidalError, simpsonError },
//     one entry per subinterval count n = 1..maxN. simpsonArea/simpsonError
//     are NaN when n is odd.
//   f — parsed function, reused for curve sampling in setup
//   trueIntegral — reference integral of f itself, shared by every row
result = {
  const math = window.math

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

  if (!f) return {rows: [], f: fallback, trueIntegral: NaN}

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

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

  // Reference integral: composite Simpson's rule with a very fine even
  // subinterval count, far finer than any n tested below.
  const trueIntegral = VM.numerical.simpsonEstimate(f, lo, hi, 4000)

  const rows = []
  let n = 1
  while (n <= maxN) {
    const h = (hi - lo) / n

    // Left Riemann sum
    let leftArea = 0
    for (let i = 0; i < n; i++) {
      const x = lo + i * h
      const y = f(x)
      if (Number.isFinite(y)) leftArea += y * h
    }

    // Midpoint Riemann sum
    let midpointArea = 0
    for (let i = 0; i < n; i++) {
      const mid = lo + (i + 0.5) * h
      const y = f(mid)
      if (Number.isFinite(y)) midpointArea += y * h
    }

    // Trapezoidal rule
    let trapezoidalArea = 0
    {
      const y0 = f(lo)
      const yn = f(hi)
      let sum = 0
      if (Number.isFinite(y0)) sum += y0 / 2
      if (Number.isFinite(yn)) sum += yn / 2
      for (let i = 1; i < n; i++) {
        const y = f(lo + i * h)
        if (Number.isFinite(y)) sum += y
      }
      trapezoidalArea = sum * h
    }

    // Simpson's rule needs an even number of subintervals.
    let simpsonArea = NaN
    if (n % 2 === 0) simpsonArea = VM.numerical.simpsonEstimate(f, lo, hi, n)

    const leftError = Math.abs(leftArea - trueIntegral)
    const midpointError = Math.abs(midpointArea - trueIntegral)
    const trapezoidalError = Math.abs(trapezoidalArea - trueIntegral)
    let simpsonError = NaN
    if (Number.isFinite(simpsonArea)) simpsonError = Math.abs(simpsonArea - trueIntegral)

    rows.push({n, leftArea, midpointArea, trapezoidalArea, simpsonArea, trueIntegral, leftError, midpointError, trapezoidalError, simpsonError})
    n++
  }

  return {rows, f, trueIntegral}
}

// setup computes the data needed by the main chart at the current n: the
// rectangles for the left Riemann sum and the midpoint Riemann sum, the
// nodes for the trapezoidal rule's linear interpolant, dense parabola
// samples for Simpson's rule (only when n is even), and the axis ranges /
// smooth f(x) curve.
setup = {
  if (result.rows.length === 0) {
    const empty = {lo: -3, hi: 3}
    return {n: 0, nodesX: [], nodesY: [], leftY: [], midX: [], midY: [], widths: [], simpsonX: [], simpsonY: [], sampleXs: [], sampleYs: [], xRange: empty, xBuffer: empty, yRange: {lo: -1, hi: 1}}
  }

  const n = Math.min(Math.max(Math.round(Number(nControl)), 1), result.rows.length)

  const a = committed.a0
  const b = committed.b0
  const lo = Math.min(a, b)
  const hi = Math.max(a, b)
  const h = (hi - lo) / n

  // Nodes for the trapezoidal rule's piecewise-linear interpolant.
  const nodesX = []
  const nodesY = []
  for (let k = 0; k <= n; k++) {
    const x = lo + k * h
    nodesX.push(x)
    nodesY.push(result.f(x))
  }

  // Left endpoint heights for the left Riemann sum's rectangles — these
  // share the same trapezoidal nodes, just dropping the last (right) one.
  const leftY = nodesY.slice(0, n)

  // Midpoints and heights for the midpoint Riemann sum's rectangles.
  const midX = []
  const midY = []
  const widths = []
  for (let i = 0; i < n; i++) {
    const mid = lo + (i + 0.5) * h
    midX.push(mid)
    midY.push(result.f(mid))
    widths.push(h)
  }

  // Dense parabola samples per pair of subintervals for Simpson's rule,
  // via the Lagrange quadratic through each triple of nodes — only defined
  // when n is even.
  const simpsonX = []
  const simpsonY = []
  if (n % 2 === 0) {
    const samplesPerPair = 40
    for (let i = 0; i < n; i += 2) {
      const panel = VM.numerical.lagrangeQuadratic(nodesX[i], nodesY[i], nodesX[i + 1], nodesY[i + 1], nodesX[i + 2], nodesY[i + 2], samplesPerPair)
      for (const x of panel.xs) simpsonX.push(x)
      for (const y of panel.ys) simpsonY.push(y)
    }
  }

  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, rectangle, and parabola 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 y of midY) if (Number.isFinite(y)) yValues.push(y)
  for (const y of simpsonY) if (Number.isFinite(y)) yValues.push(y)
  const yRange = VM.plotting.paddedRange(yValues, {emptyRange: [-1, 1], relativePadding: 0.25, minPadding: 0.5})

  return {n, nodesX, nodesY, leftY, midX, midY, widths, simpsonX, simpsonY, sampleXs, sampleYs, xRange, xBuffer, yRange}
}

// mainPlot draws f(x) and, depending on the checkboxes, the rectangles for
// the left Riemann sum and the midpoint Riemann sum, the piecewise-linear
// interpolant for the trapezoidal rule, and the piecewise-parabola
// interpolant for Simpson's rule — each shaded to show the area it
// approximates.
//
// _plotDiv is kept in a closure so Plotly reuses the same DOM node across
// reactive updates — this preserves zoom and pan when n or the checkboxes
// change.
mainPlot = {
  const Plotly = window.Plotly
  let _plotDiv = null

  return (data, showFCurve, showLeftRiemann, showMidpoint, showTrapezoidal, showSimpson, showNodes) => {
    const traces = []

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

    if (showLeftRiemann && data.midX.length > 0) {
      traces.push({
        x: data.midX, y: data.leftY, width: data.widths, type: "bar", name: "Left Riemann sum", showlegend: false,
        marker: {color: "rgba(139, 92, 246, 0.25)", line: {color: "#8b5cf6", width: 1.5}}
      })
    }

    if (showMidpoint && data.midX.length > 0) {
      traces.push({
        x: data.midX, y: data.midY, width: data.widths, type: "bar", name: "Midpoint Riemann sum", showlegend: false,
        marker: {color: VM.plotting.alpha("warn", 0.25), line: {color: chartColors.warn, width: 1.5}}
      })
    }

    if (showTrapezoidal && data.nodesX.length > 0) {
      traces.push({
        x: data.nodesX, y: data.nodesY, type: "scatter", mode: "lines", name: "Trapezoidal rule", showlegend: false,
        fill: "tozeroy", fillcolor: VM.plotting.alpha("alt", 0.2), line: {color: chartColors.alt, width: 2.5}
      })
    }

    if (showSimpson && data.simpsonX.length > 0) {
      traces.push({
        x: data.simpsonX, y: data.simpsonY, type: "scatter", mode: "lines", name: "Simpson's rule", showlegend: false,
        fill: "tozeroy", fillcolor: VM.plotting.alpha("ok", 0.2), line: {color: chartColors.ok, width: 2.5}
      })
    }

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

    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},
      barmode: "overlay",
      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
  }
}

// maxN parses the Max field's expression, falling back to 20 for invalid
// input so the n slider always gets a usable bound.
maxN = {
  const math = window.math
  const parsed = VM.expressions.makeNumber(math, maxNInput)
  if (parsed === null || !Number.isFinite(parsed)) return 20
  return Math.max(1, Math.round(parsed))
}

mainPlot(setup, showFCurve, showLeftRiemann, showMidpoint, showTrapezoidal, showSimpson, showNodes)

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

viewof nControl = {
  const nSlider = Inputs.range([1, maxN], {step: 1, value: Math.min(initialN, maxN), label: "n (subintervals)"})
  nSlider.classList.add("ojs-fill")
  nSlider.dataset.exampleField = "n"
  return nSlider
}

NoteConvergence plots

// iteratesPlot draws the left, midpoint, trapezoidal, and Simpson's area
// estimates vs n for the entire run, against a reference line at the true
// integral of f — it doesn't depend on the current step, so it's just
// plotted once per run.
iteratesPlot = {
  const leftName = "Left Riemann sum"
  const midpointName = "Midpoint Riemann sum"
  const trapezoidalName = "Trapezoidal rule"
  const simpsonName = "Simpson's rule"
  const trueName = "True integral"

  const marks = [
    Plot.ruleY([result.trueIntegral], {stroke: () => trueName, strokeDasharray: "4,3"})
  ]

  if (showLeftRiemann) {
    marks.push(
      Plot.lineY(result.rows, {x: "n", y: "leftArea", stroke: () => leftName}),
      Plot.dot(result.rows, {x: "n", y: "leftArea", fill: () => leftName, r: 3})
    )
  }

  if (showMidpoint) {
    marks.push(
      Plot.lineY(result.rows, {x: "n", y: "midpointArea", stroke: () => midpointName}),
      Plot.dot(result.rows, {x: "n", y: "midpointArea", fill: () => midpointName, r: 3})
    )
  }

  if (showTrapezoidal) {
    marks.push(
      Plot.lineY(result.rows, {x: "n", y: "trapezoidalArea", stroke: () => trapezoidalName}),
      Plot.dot(result.rows, {x: "n", y: "trapezoidalArea", fill: () => trapezoidalName, r: 3})
    )
  }

  if (showSimpson) {
    marks.push(
      Plot.lineY(result.rows, {x: "n", y: "simpsonArea", stroke: () => simpsonName}),
      Plot.dot(result.rows, {x: "n", y: "simpsonArea", fill: () => simpsonName, r: 3})
    )
  }

  const colorDomain = [leftName, midpointName, trapezoidalName, simpsonName, trueName]
  const colorRange = ["#8b5cf6", chartColors.warn, chartColors.alt, chartColors.ok, chartColors.ink]

  const chart = Plot.plot({
    width: 340, height: 300, marginLeft: 46,
    x: {label: "n"},
    y: {label: "area"},
    color: {domain: colorDomain, range: colorRange},
    marks
  })
  const legend = Plot.legend({color: {domain: colorDomain, range: colorRange}})

  const container = document.createElement("div")
  container.append(chart, legend)
  return container
}
// convergencePlot draws log₁₀|error| vs log₁₀ n for all four methods, with
// a linear regression line for each. The slope of each line estimates the
// order of convergence of that method (negated, since error shrinks as n
// grows) — for the entire run, so it's just plotted once per run.
convergencePlot = {
  const leftPts = []
  const midpointPts = []
  const trapezoidalPts = []
  const simpsonPts = []
  for (const row of result.rows) {
    if (Number.isFinite(row.leftError) && row.leftError > 0) {
      leftPts.push({x: Math.log10(row.n), y: Math.log10(row.leftError)})
    }
    if (Number.isFinite(row.midpointError) && row.midpointError > 0) {
      midpointPts.push({x: Math.log10(row.n), y: Math.log10(row.midpointError)})
    }
    if (Number.isFinite(row.trapezoidalError) && row.trapezoidalError > 0) {
      trapezoidalPts.push({x: Math.log10(row.n), y: Math.log10(row.trapezoidalError)})
    }
    if (Number.isFinite(row.simpsonError) && row.simpsonError > 0) {
      simpsonPts.push({x: Math.log10(row.n), y: Math.log10(row.simpsonError)})
    }
  }

  const leftName = "Left error"
  const midpointName = "Midpoint error"
  const trapezoidalName = "Trapezoidal error"
  const simpsonName = "Simpson's error"
  const marks = []
  const colorDomain = []
  const colorRange = []

  if (showLeftRiemann) {
    marks.push(Plot.dot(leftPts, {x: "x", y: "y", fill: () => leftName, r: 4}))
    colorDomain.push(leftName)
    colorRange.push("#8b5cf6")
  }

  if (showMidpoint) {
    marks.push(Plot.dot(midpointPts, {x: "x", y: "y", fill: () => midpointName, r: 4}))
    colorDomain.push(midpointName)
    colorRange.push(chartColors.warn)
  }

  if (showTrapezoidal) {
    marks.push(Plot.dot(trapezoidalPts, {x: "x", y: "y", fill: () => trapezoidalName, r: 4}))
    colorDomain.push(trapezoidalName)
    colorRange.push(chartColors.alt)
  }

  if (showSimpson) {
    marks.push(Plot.dot(simpsonPts, {x: "x", y: "y", fill: () => simpsonName, r: 4}))
    colorDomain.push(simpsonName)
    colorRange.push(chartColors.ok)
  }

  // legendDomain/legendRange hold only the regression lines — the legend
  // shown below the chart displays just the fitted orders, not the raw
  // point series, to keep it uncluttered.
  const legendDomain = [], legendRange = []

  const leftFit = showLeftRiemann ? VM.numerical.linearRegression(leftPts) : null
  if (leftFit) {
    const lineName = `left order ≈ ${(-leftFit.slope).toFixed(2)}`
    marks.push(Plot.line(
      [{x: leftFit.xlo, y: leftFit.slope * leftFit.xlo + leftFit.intercept}, {x: leftFit.xhi, y: leftFit.slope * leftFit.xhi + leftFit.intercept}],
      {x: "x", y: "y", stroke: () => lineName, strokeDasharray: "6,4"}
    ))
    colorDomain.push(lineName)
    colorRange.push("#a78bfa")
    legendDomain.push(lineName)
    legendRange.push("#a78bfa")
  }

  const midpointFit = showMidpoint ? VM.numerical.linearRegression(midpointPts) : null
  if (midpointFit) {
    const lineName = `midpoint order ≈ ${(-midpointFit.slope).toFixed(2)}`
    marks.push(Plot.line(
      [{x: midpointFit.xlo, y: midpointFit.slope * midpointFit.xlo + midpointFit.intercept}, {x: midpointFit.xhi, y: midpointFit.slope * midpointFit.xhi + midpointFit.intercept}],
      {x: "x", y: "y", stroke: () => lineName, strokeDasharray: "6,4"}
    ))
    colorDomain.push(lineName)
    colorRange.push("#fbbf24")
    legendDomain.push(lineName)
    legendRange.push("#fbbf24")
  }

  const trapezoidalFit = showTrapezoidal ? VM.numerical.linearRegression(trapezoidalPts) : null
  if (trapezoidalFit) {
    const lineName = `trapezoidal order ≈ ${(-trapezoidalFit.slope).toFixed(2)}`
    marks.push(Plot.line(
      [{x: trapezoidalFit.xlo, y: trapezoidalFit.slope * trapezoidalFit.xlo + trapezoidalFit.intercept}, {x: trapezoidalFit.xhi, y: trapezoidalFit.slope * trapezoidalFit.xhi + trapezoidalFit.intercept}],
      {x: "x", y: "y", stroke: () => lineName, strokeDasharray: "6,4"}
    ))
    colorDomain.push(lineName)
    colorRange.push("#f87171")
    legendDomain.push(lineName)
    legendRange.push("#f87171")
  }

  const simpsonFit = showSimpson ? VM.numerical.linearRegression(simpsonPts) : null
  if (simpsonFit) {
    const lineName = `simpson order ≈ ${(-simpsonFit.slope).toFixed(2)}`
    marks.push(Plot.line(
      [{x: simpsonFit.xlo, y: simpsonFit.slope * simpsonFit.xlo + simpsonFit.intercept}, {x: simpsonFit.xhi, y: simpsonFit.slope * simpsonFit.xhi + simpsonFit.intercept}],
      {x: "x", y: "y", stroke: () => lineName, strokeDasharray: "6,4"}
    ))
    colorDomain.push(lineName)
    colorRange.push("#4ade80")
    legendDomain.push(lineName)
    legendRange.push("#4ade80")
  }

  const chart = Plot.plot({
    width: 340, height: 300, marginLeft: 46,
    x: {label: "log₁₀ n"},
    y: {label: "log₁₀ |error|"},
    color: {domain: colorDomain, range: colorRange},
    marks
  })

  const container = document.createElement("div")
  container.append(chart)
  // Plot.legend throws on an empty domain/range — there's no legend to
  // show when none of the checked methods produced a regression line
  // (e.g. too few finite error points to fit).
  if (legendDomain.length > 0) {
    const legend = Plot.legend({color: {domain: legendDomain, range: legendRange}})
    container.append(legend)
  }
  return container
}

NoteArea table

// Area table — shows the entire run, not just the current n, with all
// four method estimates and their errors against the true integral.
areaTable = {
  if (result.rows.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 (const row of result.rows) {
    let simpsonCell = "NaN"
    if (Number.isFinite(row.simpsonArea)) simpsonCell = row.simpsonArea.toPrecision(10)

    let simpsonErrorCell = "NaN"
    if (Number.isFinite(row.simpsonError)) simpsonErrorCell = row.simpsonError.toExponential(4)

    formattedRows.push([
      row.n,
      row.leftArea.toPrecision(10),
      row.midpointArea.toPrecision(10),
      row.trapezoidalArea.toPrecision(10),
      simpsonCell,
      row.trueIntegral.toPrecision(10),
      row.leftError.toExponential(4),
      row.midpointError.toExponential(4),
      row.trapezoidalError.toExponential(4),
      simpsonErrorCell
    ])
  }

  return VM.ui.renderTable({
    html,
    headers: ["n", tex`L_n`, tex`M_n`, tex`T_n`, tex`S_n`, tex`I`, tex`|L_n - I|`, tex`|M_n - I|`, tex`|T_n - I|`, tex`|S_n - I|`],
    csvHeaders: ["n", "L_n (left)", "M_n (midpoint)", "T_n (trapezoidal)", "S_n (simpson)", "I (true integral)", "|L_n - I|", "|M_n - I|", "|T_n - I|", "|S_n - I|"],
    filename: "integration-methods.csv",
    rows: formattedRows
  })
}

Approximates \(\int_a^b f(x)\, dx\) by replacing \(f\) with a simple piecewise interpolant on \(n\) equally spaced subintervals of width \(h=(b-a)/n\), at nodes \(x_i=a+ih\), and integrating that interpolant exactly:

\[ L_n = h\sum_{i=0}^{n-1} f(x_i), \qquad M_n = h\sum_{i=0}^{n-1} f\!\left(x_i+\tfrac{h}{2}\right), \]

\[ T_n = h\left(\tfrac{1}{2}f(x_0)+f(x_1)+\cdots+f(x_{n-1})+\tfrac{1}{2}f(x_n)\right), \]

\[ S_n = \frac{h}{3}\Big(f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+\cdots+4f(x_{n-1})+f(x_n)\Big)\quad(n\text{ even}), \]

a left or midpoint rectangle sum, a trapezoid sum, or (for Simpson’s rule) a quadratic through each pair of subintervals.