Linear vs. Cubic Spline Interpolation

// VM is the shared utility library, loaded globally via _includes/head-scripts.html (js/**), plus this page's own cubic-spline.js (loaded via the include-in-header above). This page uses VM.cubicSpline, VM.expressions.makeFunction, VM.expressions.makeNumber, VM.plotting.paddedRange, VM.ui.renderTable.
VM = window.VM

Given a function \(f\) on \([a, b]\), sample it at \(n+1\) equally spaced nodes \(x_0 < x_1 < \dots < x_n\) and connect the samples \((x_i, f(x_i))\) with simple pieces. A linear interpolant joins consecutive samples with straight segments; a cubic spline interpolant joins them with cubic polynomials chosen so the curve and its first two derivatives match up at each node.

As \(n\) grows, both interpolants should converge to \(f\) — but at different rates. This page tracks that convergence through the arc length of each interpolant,

\[ L = \int_a^b \sqrt{1 + g'(x)^2}\, dx, \]

comparing \(L_{\text{linear}}(n)\) and \(L_{\text{cubic}}(n)\) against the true arc length of \(f\) as \(n\) increases.

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") ?? "2*pi", placeholder: "Example: 2*pi, 6.28"})
  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 interpolant 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 interpolants = {
  const options = ["Linear interpolant", "Cubic spline interpolant"]
  const checkboxInput = Inputs.checkbox(options, {value: options})
  checkboxInput.classList.add("ojs-fill")
  return checkboxInput
}

showLinear = interpolants.includes("Linear interpolant")

showCubic = interpolants.includes("Cubic spline interpolant")

// 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 arc-length
// history used by the convergence plots and tables.
//   rows — array of { n, linearLength, cubicLength, trueArcLength, linearError, cubicError },
//     one entry per subinterval count n = 1..maxN
//   f — parsed function, reused for curve sampling in setup
//   trueArcLength — reference arc length 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, trueArcLength: NaN}

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

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

  // Reference arc length of f: chord-length sum over a very fine uniform
  // grid, far finer than any interpolant tested below.
  const refCount = 4000
  let trueArcLength = 0
  let xPrev = lo
  let yPrev = f(xPrev)
  for (let k = 1; k <= refCount; k++) {
    const x = lo + (hi - lo) * k / refCount
    const y = f(x)
    if (Number.isFinite(yPrev) && Number.isFinite(y)) {
      const dx = x - xPrev
      const dy = y - yPrev
      trueArcLength += Math.sqrt(dx * dx + dy * dy)
    }
    xPrev = x
    yPrev = y
  }

  const rows = []
  let n = 1
  while (n <= maxN) {
    // n+1 equally spaced interpolation nodes over [lo, hi]
    const nodesX = []
    const nodesY = []
    for (let k = 0; k <= n; k++) {
      const x = lo + (hi - lo) * k / n
      nodesX.push(x)
      nodesY.push(f(x))
    }

    // Linear interpolant arc length: exact sum of segment lengths.
    let linearLength = 0
    for (let k = 0; k < n; k++) {
      const dx = nodesX[k + 1] - nodesX[k]
      const dy = nodesY[k + 1] - nodesY[k]
      if (Number.isFinite(dy)) linearLength += Math.sqrt(dx * dx + dy * dy)
    }

    // Cubic spline arc length: approximated by sampling the spline densely
    // and summing chord lengths along the sampled curve.
    const spline = VM.cubicSpline(nodesX, nodesY)
    let cubicLength = NaN
    if (spline) {
      const sampleCount = n * 20
      let sxPrev = nodesX[0]
      let syPrev = spline(sxPrev)
      cubicLength = 0
      for (let k = 1; k <= sampleCount; k++) {
        const x = lo + (hi - lo) * k / sampleCount
        const y = spline(x)
        if (Number.isFinite(syPrev) && Number.isFinite(y)) {
          const dx = x - sxPrev
          const dy = y - syPrev
          cubicLength += Math.sqrt(dx * dx + dy * dy)
        }
        sxPrev = x
        syPrev = y
      }
    }

    const linearError = Math.abs(linearLength - trueArcLength)
    const cubicError = Number.isFinite(cubicLength) ? Math.abs(cubicLength - trueArcLength) : NaN

    rows.push({n, linearLength, cubicLength, trueArcLength, linearError, cubicError})
    n++
  }

  return {rows, f, trueArcLength}
}

// setup computes the data needed by the main chart at the current n: the
// interpolation nodes, the linear polyline, a dense sample of the cubic
// spline, 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: [], linearX: [], linearY: [], cubicX: [], cubicY: [], 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)

  // Interpolation nodes for the current n.
  const nodesX = []
  const nodesY = []
  for (let k = 0; k <= n; k++) {
    const x = lo + (hi - lo) * k / n
    nodesX.push(x)
    nodesY.push(result.f(x))
  }

  // Linear interpolant is just the node polyline.
  const linearX = nodesX.slice()
  const linearY = nodesY.slice()

  // Cubic spline: dense samples for a smooth curve.
  const spline = VM.cubicSpline(nodesX, nodesY)
  const cubicX = []
  const cubicY = []
  if (spline) {
    const sampleCount = Math.max(200, n * 20)
    for (let k = 0; k <= sampleCount; k++) {
      const x = lo + (hi - lo) * k / sampleCount
      cubicX.push(x)
      cubicY.push(spline(x))
    }
  }

  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 cubic-spline values.
  const yValues = []
  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 cubicY) 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, linearX, linearY, cubicX, cubicY, sampleXs, sampleYs, xRange, xBuffer, yRange}
}

// mainPlot draws f(x) and, depending on the checkboxes, the linear and/or
// cubic spline interpolant through the current n's nodes.
//
// _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, showLinear, showCubic) => {
    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 (showLinear && data.linearX.length > 0) {
      traces.push({x: data.linearX, y: data.linearY, type: "scatter", mode: "lines", name: "Linear interpolant", showlegend: true, line: {color: "#dc2626", width: 2.5}})
    }
    if (showCubic && data.cubicX.length > 0) {
      traces.push({x: data.cubicX, y: data.cubicY, type: "scatter", mode: "lines", name: "Cubic spline interpolant", showlegend: true, line: {color: "#16a34a", width: 2.5}})
    }
    if (data.nodesX.length > 0 && (showLinear || showCubic)) {
      traces.push({x: data.nodesX, y: data.nodesY, type: "scatter", mode: "markers", name: "Interpolation nodes", showlegend: true, marker: {color: "#111827", size: 8, 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},
      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, showLinear, showCubic)

NoteConvergence plots

// iteratesPlot draws the linear and cubic spline arc lengths vs n for the
// entire run, against a reference line at the true arc length of f — it
// doesn't depend on the current step, so it's just plotted once per run.
iteratesPlot = {
  const linearName = "Linear interpolant"
  const cubicName = "Cubic spline interpolant"
  const trueName = "True arc length"

  const marks = [
    Plot.ruleY([result.trueArcLength], {stroke: () => trueName, strokeDasharray: "4,3"}),
    Plot.lineY(result.rows, {x: "n", y: "linearLength", stroke: () => linearName}),
    Plot.dot(result.rows, {x: "n", y: "linearLength", fill: () => linearName, r: 3}),
    Plot.lineY(result.rows, {x: "n", y: "cubicLength", stroke: () => cubicName}),
    Plot.dot(result.rows, {x: "n", y: "cubicLength", fill: () => cubicName, r: 3})
  ]
  const colorDomain = [linearName, cubicName, trueName]
  const colorRange = ["#dc2626", "#16a34a", "#111827"]

  const chart = Plot.plot({
    width: 340, height: 300, marginLeft: 46,
    x: {label: "n"},
    y: {label: "arc length"},
    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 both interpolants, with
// a linear regression line for each. The slope of each line estimates the
// order of convergence of that interpolant's arc length (negated, since
// error shrinks as n grows) — for the entire run, so it's just plotted once
// per run.
convergencePlot = {
  const regression = points => {
    if (points.length < 2) return null
    const pointCount = points.length
    let sx = 0, sy = 0, sxy = 0, sxx = 0
    for (const p of points) {
      sx += p.x
      sy += p.y
      sxy += p.x * p.y
      sxx += p.x * p.x
    }
    const denom = pointCount * sxx - sx * sx
    if (denom === 0) return null
    const slope = (pointCount * sxy - sx * sy) / denom
    const intercept = (sy - slope * sx) / pointCount
    return {slope, intercept, xlo: points[0].x, xhi: points[pointCount - 1].x}
  }

  const linearPts = []
  const cubicPts = []
  for (const row of result.rows) {
    if (Number.isFinite(row.linearError) && row.linearError > 0) {
      linearPts.push({x: Math.log10(row.n), y: Math.log10(row.linearError)})
    }
    if (Number.isFinite(row.cubicError) && row.cubicError > 0) {
      cubicPts.push({x: Math.log10(row.n), y: Math.log10(row.cubicError)})
    }
  }

  const linearName = "Linear interpolant error"
  const cubicName = "Cubic spline error"
  const marks = [
    Plot.dot(linearPts, {x: "x", y: "y", fill: () => linearName, r: 4}),
    Plot.dot(cubicPts, {x: "x", y: "y", fill: () => cubicName, r: 4})
  ]
  const colorDomain = [linearName, cubicName]
  const colorRange = ["#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 linearFit = regression(linearPts)
  if (linearFit) {
    const lineName = `linear order ≈ ${(-linearFit.slope).toFixed(2)}`
    marks.push(Plot.line(
      [{x: linearFit.xlo, y: linearFit.slope * linearFit.xlo + linearFit.intercept}, {x: linearFit.xhi, y: linearFit.slope * linearFit.xhi + linearFit.intercept}],
      {x: "x", y: "y", stroke: () => lineName, strokeDasharray: "6,4"}
    ))
    colorDomain.push(lineName)
    colorRange.push("#f87171")
    legendDomain.push(lineName)
    legendRange.push("#f87171")
  }

  const cubicFit = regression(cubicPts)
  if (cubicFit) {
    const lineName = `cubic order ≈ ${(-cubicFit.slope).toFixed(2)}`
    marks.push(Plot.line(
      [{x: cubicFit.xlo, y: cubicFit.slope * cubicFit.xlo + cubicFit.intercept}, {x: cubicFit.xhi, y: cubicFit.slope * cubicFit.xhi + cubicFit.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
}

NoteArc length table

// Arc length table — shows the entire run, not just the current n, with
// both arc lengths and both errors against the true arc length of f.
arcLengthTable = {
  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 cubicCell = "NaN"
    if (Number.isFinite(row.cubicLength)) cubicCell = row.cubicLength.toPrecision(10)

    let cubicErrorCell = "NaN"
    if (Number.isFinite(row.cubicError)) cubicErrorCell = row.cubicError.toExponential(4)

    formattedRows.push([
      row.n,
      row.linearLength.toPrecision(10),
      cubicCell,
      row.trueArcLength.toPrecision(10),
      row.linearError.toExponential(4),
      cubicErrorCell
    ])
  }

  return VM.ui.renderTable({
    html,
    headers: ["n", tex`L_{\text{linear}}`, tex`L_{\text{cubic}}`, tex`L_{\text{true}}`, tex`|L_{\text{linear}} - L_{\text{true}}|`, tex`|L_{\text{cubic}} - L_{\text{true}}|`],
    csvHeaders: ["n", "L_linear", "L_cubic", "L_true", "|L_linear - L_true|", "|L_cubic - L_true|"],
    filename: "linear-cubic-arc-lengths.csv",
    rows: formattedRows
  })
}