Explicit Methods for ODEs

Compare Euler, midpoint and RK4 on the same initial value problem.

// VM is the shared utility library, loaded on every page by the mathviz extension. This page uses VM.numerical.eulerSolve, VM.numerical.linearRegression, VM.expressions.makeFunction2, VM.expressions.makeNumber, VM.plotting.paddedRange, VM.ui.renderTable, VM.numerical.rk4Solve, VM.ui.applyExampleParams.
VM = window.VM

urlParams = new URLSearchParams(window.location.search)

methodOptions = ["Reference solution", "Euler's method", "Improved Euler", "RK4"]

// Matches the line/marker colors mainPlot uses for each method's 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)

methodColors = new Map([
  ["Reference solution", chartColors.fn],
  ["Euler's method", chartColors.alt],
  ["Improved Euler", chartColors.warn],
  ["RK4", chartColors.ok]
])

methodSlugToLabel = new Map([
  ["reference", "Reference solution"],
  ["reference-solution", "Reference solution"],
  ["ref", "Reference solution"],
  ["true", "Reference solution"],
  ["true-solution", "Reference solution"],
  ["exact", "Reference solution"],

  ["euler", "Euler's method"],
  ["eulers", "Euler's method"],
  ["euler-method", "Euler's method"],
  ["eulers-method", "Euler's method"],
  ["euler_method", "Euler's method"],

  ["improved", "Improved Euler"],
  ["improved-euler", "Improved Euler"],
  ["improvedeuler", "Improved Euler"],
  ["improved_euler", "Improved Euler"],

  ["rk4", "RK4"],
  ["runge-kutta", "RK4"],
  ["runge-kutta-4", "RK4"],
  ["rk-4", "RK4"]
])

methodLabelToSlug = new Map([
  ["Reference solution", "reference"],
  ["Euler's method", "euler"],
  ["Improved Euler", "improved"],
  ["RK4", "rk4"]
])

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) : 3
}

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

// The page's examples and method live in method.js (the fit is the
// library's VM.numerical.orderFromErrorsVsN), so
// method.test.js can run every example and check what the convergence plot
// reports under npm test.
import {examples, methodErrors, improvedEulerSolve, endpointErrors} from "./method.js"

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

// 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 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 y0 = {
  const params = new URLSearchParams(window.location.search)
  const y0Input = Inputs.text({label: "Initial condition y(0)", value: params.get("y0") ?? "1", placeholder: "Example: 0.5, 1/2"})
  y0Input.dataset.exampleField = "y0"
  return y0Input
}

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

// 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 : "Custom inputs", 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"]',
  y0: '[data-example-field="y0"]',
  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)
  })
}

// Keeps the dropdown naming the example the fields above hold ("Custom
// inputs" once none does). Reads every example field's value so it re-runs
// whenever one changes; it only sets the dropdown, never a field.
syncExampleFromFields = {
  fText; y0; bEnd; maxNInput; nControl
  VM.ui.syncExampleSelect(viewof exampleSelect, examples, exampleFieldSelectors)
}

showRefSolution = methods.includes("Reference solution")

showEuler = methods.includes("Euler's method")

showImprovedEuler = methods.includes("Improved Euler")

showRK4 = methods.includes("RK4")

// 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/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)
  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, y0: y0Val, bEnd: bVal}
}

// result: for the current committed function/initial condition/endpoint,
// builds the endpoint-value history used by the convergence plots and
// table.
//   rows — array of { n, eulerValue, improvedEulerValue, rk4Value, trueValue,
//     eulerError, improvedEulerError, rk4Error }, one entry per subinterval
//     count n = 1..maxN — the computed y(b) for each method against a
//     reference value.
//   f — parsed right-hand side, reused for the trajectories in setup
//   trueValue — reference y(b), shared by every row
result = {
  const math = window.math

  const f = VM.expressions.makeFunction2(math, committed.fText)
  const fallback = (t, y) => NaN

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

  const t0 = 0
  const y0 = committed.y0
  const b = committed.bEnd
  if (!Number.isFinite(y0) || !Number.isFinite(b) || b === t0) return {rows: [], f, trueValue: NaN}

  const {rows, trueValue} = endpointErrors(f, y0, b, maxN)
  return {rows, f, trueValue}
}

// setup computes the data needed by the main chart at the current n: each
// method's step-by-step trajectory, and a fine reference trajectory for
// the smooth "true" solution curve.
setup = {
  if (result.rows.length === 0) {
    const empty = {lo: -3, hi: 3}
    return {n: 0, eulerTs: [], eulerYs: [], improvedEulerTs: [], improvedEulerYs: [], rk4Ts: [], rk4Ys: [], refTs: [], refYs: [], tRange: empty, yRange: {lo: -1, hi: 1}}
  }

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

  const t0 = 0
  const y0 = committed.y0
  const b = committed.bEnd

  const euler = VM.numerical.eulerSolve(result.f, t0, y0, b, n)
  const improvedEuler = improvedEulerSolve(result.f, t0, y0, b, n)
  const rk4 = VM.numerical.rk4Solve(result.f, t0, y0, b, n)

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

  // The y-range fits the reference solution only, not the methods'
  // trajectories, so moving the n slider never rescales the axis — a method
  // that is off at small n visibly runs off the chart instead. Cap at 1e8
  // so a blowing-up reference doesn't wreck the range.
  const yValues = []
  for (const y of reference.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})

  return {
    n,
    eulerTs: euler.ts, eulerYs: euler.ys,
    improvedEulerTs: improvedEuler.ts, improvedEulerYs: improvedEuler.ys,
    rk4Ts: rk4.ts, rk4Ys: rk4.ys,
    refTs: reference.ts, refYs: reference.ys,
    tRange, yRange
  }
}

// mainPlot draws the fine reference solution curve and, depending on the
// checkboxes, the piecewise-linear trajectory computed by each method at
// the current n.
//
// _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, showRefSolution, showEuler, showImprovedEuler, showRK4) => {
    const traces = []

    if (showRefSolution) {
      traces.push({x: data.refTs, y: data.refYs, type: "scatter", mode: "lines", name: "Reference solution", showlegend: false, line: {color: chartColors.fn, width: 3}})
    }

    if (showEuler && data.eulerTs.length > 0) {
      traces.push({
        x: data.eulerTs, y: data.eulerYs, type: "scatter", mode: "lines+markers", name: "Euler's method", showlegend: false,
        line: {color: chartColors.alt, width: 2}, marker: {color: chartColors.alt, size: 6}
      })
    }
    if (showImprovedEuler && data.improvedEulerTs.length > 0) {
      traces.push({
        x: data.improvedEulerTs, y: data.improvedEulerYs, type: "scatter", mode: "lines+markers", name: "Improved Euler", showlegend: false,
        line: {color: chartColors.warn, width: 2}, marker: {color: chartColors.warn, size: 6}
      })
    }
    if (showRK4 && data.rk4Ts.length > 0) {
      traces.push({
        x: data.rk4Ts, y: data.rk4Ys, type: "scatter", mode: "lines+markers", name: "RK4", showlegend: false,
        line: {color: chartColors.ok, width: 2}, marker: {color: chartColors.ok, size: 6}
      })
    }

    const layout = {
      xaxis: {title: "t", range: [data.tRange.lo, data.tRange.hi], zeroline: true},
      yaxis: {title: "y", range: [data.yRange.lo, data.yRange.hi], zeroline: true},
      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, showRefSolution, showEuler, showImprovedEuler, showRK4)

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
}

Given \(y' = f(t, y)\), \(y(0) = y_0\), an explicit method computes each \(y_{i+1} \approx y(t_{i+1})\) directly from already-known values, using \(h\) for the step size and \(t_i = t_0 + ih\). Compare Euler’s method (one slope evaluation per step):

\[ y_{i+1} = y_i + h\,f(t_i, y_i), \]

Improved Euler (averages the slope at both endpoints of the step):

\[ \begin{aligned} k_1 &= f(t_i, y_i), \\ k_2 &= f(t_{i+1},\, y_i + h k_1), \\ y_{i+1} &= y_i + \frac{h}{2}(k_1 + k_2), \end{aligned} \]

and RK4 (averages four slope estimates):

\[ \begin{aligned} k_1 &= f(t_i, y_i), \\ k_2 &= f\!\left(t_i + \tfrac{h}{2},\, y_i + \tfrac{h}{2}k_1\right), \\ k_3 &= f\!\left(t_i + \tfrac{h}{2},\, y_i + \tfrac{h}{2}k_2\right), \\ k_4 &= f(t_{i+1},\, y_i + h k_3), \\ y_{i+1} &= y_i + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4). \end{aligned} \]

More evaluations per step buys higher accuracy.

When f depends on t alone

If \(f(t, y) = g(t)\) has no \(y\) in it, solving the ODE is integration: \(y(b) = y_0 + \int_0^b g(t)\,dt\). Each method then becomes a quadrature rule from Numerical Integration, applied on the grid \(t_i\):

What each method does when \(f(t, y) = g(t)\).

Method	Adds each step	Quadrature rule	Error
Euler	\(h\,g(t_i)\)	Left Riemann sum	\(O(h)\)
Improved Euler	\(\frac{h}{2}\big(g(t_i) + g(t_{i+1})\big)\)	Trapezoidal rule	\(O(h^2)\)
RK4	\(\frac{h}{6}\big(g(t_i) + 4g(t_i + \tfrac{h}{2}) + g(t_{i+1})\big)\)	Simpson’s rule	\(O(h^4)\)

For RK4, \(k_2 = k_3 = g(t_i + \tfrac{h}{2})\), so \(n\) steps of RK4 are Simpson’s rule on \(2n\) subintervals.

For example, the equation \(y' = 4/(1 + t^2)\), \(y(0) = 0\) (try it) has \(y(1) = \pi\), the same number as \(\int_0^1 \frac{4}{1+t^2}\,dt\).

NoteConvergence plots

// iteratesPlot draws each method's computed y(b) vs n for the entire run,
// against a reference line at the reference value of y(b) — it doesn't
// depend on the current step, so it's just plotted once per run.
iteratesPlot = {
  const eulerName = "Euler's method"
  const improvedEulerName = "Improved Euler"
  const rk4Name = "RK4"
  const trueName = "Reference value"

  const marks = [
    Plot.ruleY([result.trueValue], {stroke: () => trueName, strokeDasharray: "4,3"}),
    Plot.lineY(result.rows, {x: "n", y: "eulerValue", stroke: () => eulerName}),
    Plot.dot(result.rows, {x: "n", y: "eulerValue", fill: () => eulerName, r: 3}),
    Plot.lineY(result.rows, {x: "n", y: "improvedEulerValue", stroke: () => improvedEulerName}),
    Plot.dot(result.rows, {x: "n", y: "improvedEulerValue", fill: () => improvedEulerName, r: 3}),
    Plot.lineY(result.rows, {x: "n", y: "rk4Value", stroke: () => rk4Name}),
    Plot.dot(result.rows, {x: "n", y: "rk4Value", fill: () => rk4Name, r: 3})
  ]
  const colorDomain = [eulerName, improvedEulerName, rk4Name, trueName]
  const colorRange = [chartColors.alt, chartColors.warn, chartColors.ok, chartColors.ink]

  const chart = Plot.plot({
    width: 340, height: 300, marginLeft: 46,
    x: {label: "n"},
    y: {label: "y(b)"},
    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 three 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.
//
// The fit itself (which points count, and when a method is reported
// exact) is VM.numerical.orderFromErrorsVsN.
convergencePlot = {
  const methodSpecs = [
    {name: "Euler error", color: chartColors.alt, fitColor: "#f87171"},
    {name: "Improved Euler error", color: chartColors.warn, fitColor: "#fbbf24"},
    {name: "RK4 error", color: chartColors.ok, fitColor: "#4ade80"}
  ]

  const marks = []
  const colorDomain = []
  const colorRange = []
  // legendDomain/legendRange hold only the regression lines (and "exact"
  // entries) — the legend shown below the chart displays just the fitted
  // orders, not the raw point series, to keep it uncluttered.
  const legendDomain = []
  const legendRange = []

  for (let i = 0; i < methodSpecs.length; i++) {
    const spec = methodSpecs[i]
    const method = methodErrors[i]
    const {points, fit, order, exact, exactFromN} = VM.numerical.orderFromErrorsVsN(result.rows, result.trueValue, method.errorKey)

    marks.push(Plot.dot(points, {x: "x", y: "y", fill: () => spec.name, r: 4}))
    colorDomain.push(spec.name)
    colorRange.push(spec.color)

    let legendName = null
    if (exact) {
      legendName = `${method.name}: exact (error at rounding level)`
    } else if (fit) {
      legendName = `${method.name} order ≈ ${order.toFixed(2)}`
      marks.push(Plot.line(
        [{x: fit.xlo, y: fit.slope * fit.xlo + fit.intercept}, {x: fit.xhi, y: fit.slope * fit.xhi + fit.intercept}],
        {x: "x", y: "y", stroke: () => legendName, strokeDasharray: "6,4"}
      ))
    } else if (exactFromN !== null) {
      legendName = `${method.name}: exact from n = ${exactFromN}`
    }
    if (legendName === null) continue
    colorDomain.push(legendName)
    colorRange.push(spec.fitColor)
    legendDomain.push(legendName)
    legendRange.push(spec.fitColor)
  }

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

NoteValue table

// Value table — shows the entire run, not just the current n, with all
// three method estimates of y(b) and their errors against the reference
// value.
valueTable = {
  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) {
    formattedRows.push([
      row.n,
      row.eulerValue.toPrecision(10),
      row.improvedEulerValue.toPrecision(10),
      row.rk4Value.toPrecision(10),
      row.trueValue.toPrecision(10),
      row.eulerError.toExponential(4),
      row.improvedEulerError.toExponential(4),
      row.rk4Error.toExponential(4)
    ])
  }

  return VM.ui.renderTable({
    html,
    headers: ["n", tex`y_{\text{Euler}}(b)`, tex`y_{\text{Imp. Euler}}(b)`, tex`y_{\text{RK4}}(b)`, tex`y_{\text{ref}}(b)`, tex`|y_{\text{Euler}} - y_{\text{ref}}|`, tex`|y_{\text{Imp. Euler}} - y_{\text{ref}}|`, tex`|y_{\text{RK4}} - y_{\text{ref}}|`],
    csvHeaders: ["n", "y_Euler(b)", "y_ImprovedEuler(b)", "y_RK4(b)", "y_ref(b)", "|y_Euler - y_ref|", "|y_ImprovedEuler - y_ref|", "|y_RK4 - y_ref|"],
    filename: "explicit-methods-values.csv",
    rows: formattedRows
  })
}