Numerical Integration

Computational Math

Integration

Riemann Sums

Trapezoidal Rule

Simpson's Rule

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

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

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"})
  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"})
  return b0Input
}

// 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 checkboxes 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 nControl = {
  const nSlider = Inputs.range([1, maxN], {step: 1, value: Math.min(4, maxN), label: "n (subintervals)"})
  nSlider.classList.add("ojs-fill")
  return nSlider
}

viewof maxNInput = {
  const input = Inputs.text({value: "20", label: "Max", placeholder: "Example: 20, 5*4"})
  input.classList.add("ojs-auto")
  return input
}

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

viewof methods = {
  const options = ["Left Riemann sum", "Midpoint Riemann sum", "Trapezoidal rule", "Simpson's rule"]
  const checkboxInput = Inputs.checkbox(options, {value: options})
  checkboxInput.classList.add("ojs-fill")
  return checkboxInput
}

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

// 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, showLeftRiemann, showMidpoint, showTrapezoidal, showSimpson) => {
    const traces = [
      // The function curve
      {x: data.sampleXs, y: data.sampleYs, type: "scatter", mode: "lines", name: "f(x)", showlegend: true, line: {color: "#2563eb", 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: true,
        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: true,
        marker: {color: "rgba(245, 158, 11, 0.25)", line: {color: "#f59e0b", 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: true,
        fill: "tozeroy", fillcolor: "rgba(220, 38, 38, 0.2)", line: {color: "#dc2626", 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: true,
        fill: "tozeroy", fillcolor: "rgba(22, 163, 74, 0.2)", line: {color: "#16a34a", width: 2.5}
      })
    }

    if (data.nodesX.length > 0 && (showTrapezoidal || showSimpson)) {
      traces.push({x: data.nodesX, y: data.nodesY, type: "scatter", mode: "markers", name: "Interpolation nodes", showlegend: true, marker: {color: "#111827", size: 7, symbol: "circle", line: {color: "white", width: 1}}})
    }

    const layout = {
      margin: {l: 0, r: 0, t: 0, b: 0},
      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",
      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
  }
}

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

Convergence 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"}),
    Plot.lineY(result.rows, {x: "n", y: "leftArea", stroke: () => leftName}),
    Plot.dot(result.rows, {x: "n", y: "leftArea", fill: () => leftName, r: 3}),
    Plot.lineY(result.rows, {x: "n", y: "midpointArea", stroke: () => midpointName}),
    Plot.dot(result.rows, {x: "n", y: "midpointArea", fill: () => midpointName, r: 3}),
    Plot.lineY(result.rows, {x: "n", y: "trapezoidalArea", stroke: () => trapezoidalName}),
    Plot.dot(result.rows, {x: "n", y: "trapezoidalArea", fill: () => trapezoidalName, r: 3}),
    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", "#f59e0b", "#dc2626", "#16a34a", "#111827"]

  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 = [
    Plot.dot(leftPts, {x: "x", y: "y", fill: () => leftName, r: 4}),
    Plot.dot(midpointPts, {x: "x", y: "y", fill: () => midpointName, r: 4}),
    Plot.dot(trapezoidalPts, {x: "x", y: "y", fill: () => trapezoidalName, r: 4}),
    Plot.dot(simpsonPts, {x: "x", y: "y", fill: () => simpsonName, r: 4})
  ]
  const colorDomain = [leftName, midpointName, trapezoidalName, simpsonName]
  const colorRange = ["#8b5cf6", "#f59e0b", "#dc2626", "#16a34a"]
  // 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 = VM.numerical.linearRegression(leftPts)
  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 = VM.numerical.linearRegression(midpointPts)
  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 = VM.numerical.linearRegression(trapezoidalPts)
  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 = VM.numerical.linearRegression(simpsonPts)
  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 legend = Plot.legend({color: {domain: legendDomain, range: legendRange}})

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

Area 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 and integrating that interpolant exactly — a left or midpoint rectangle, a trapezoid, or (for Simpson’s rule) a quadratic through each pair of subintervals.