Implicit Methods for ODEs

// VM is the shared utility library, loaded globally via _includes/head-scripts.html (js/**). 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 = window.VM

Given an initial value problem \(y' = f(t, y)\), \(y(0) = y_0\), Euler’s method (forward/explicit Euler) uses the slope at the start of the step:

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

Implicit Euler (backward Euler) uses the slope at the end of the step instead:

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

Because \(y_{i+1}\) appears on both sides, each step requires solving a root-finding problem \(g(y_{i+1}) = y_{i+1} - y_i - h\, f(t_{i+1}, y_{i+1}) = 0\). This page solves it with a few steps of Newton’s method, starting from the explicit Euler step as the initial guess and using a numerical estimate of \(\partial f/\partial y\) in place of \(g'\).

This extra work at each step buys stability: implicit Euler stays well-behaved even with a large step size \(h\) on equations where explicit Euler oscillates or blows up. The default example below, \(y' = -5y\), is chosen to make this difference visible — try dragging \(n\) down to see explicit Euler destabilize while implicit Euler stays smooth.

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

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: 1, 1/2"})
  return y0Input
}

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

// 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
// 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 = ["Euler's method", "Implicit Euler"]
  const checkboxInput = Inputs.checkbox(options, {value: options})
  checkboxInput.classList.add("ojs-fill")
  return checkboxInput
}

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

showImplicitEuler = methods.includes("Implicit Euler")

// 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)
  history.replaceState(null, "", `${window.location.pathname}?${params}`)

  return {fText: fVal, y0: y0Val, bEnd: bVal}
}

// partialY(f, t, y) numerically estimates ∂f/∂y(t, y) via a central
// difference, used by the Newton solve inside implicitEulerSolve.
partialY = (f, t, y) => {
  const eps = 1e-6 * Math.max(1, Math.abs(y))
  return (f(t, y + eps) - f(t, y - eps)) / (2 * eps)
}

// implicitEulerSolve(f, t0, y0, tEnd, n) advances y' = f(t, y), y(t0) = y0
// with n steps of implicit (backward) Euler: y_{i+1} = y_i + h f(t_{i+1}, y_{i+1}).
// Each step solves this root-finding problem for y_{i+1} with a few
// iterations of Newton's method, starting from the explicit Euler step as
// the initial guess. Returns the full trajectory {ts, ys}.
implicitEulerSolve = (f, t0, y0, tEnd, n) => {
  const h = (tEnd - t0) / n
  const ts = [t0]
  const ys = [y0]
  let t = t0, y = y0
  for (let i = 0; i < n; i++) {
    const tNext = t0 + (i + 1) * h

    // Newton's method for g(z) = z - y - h f(tNext, z) = 0.
    let z = y + h * f(t, y)   // explicit Euler step as the initial guess
    for (let k = 0; k < 8; k++) {
      const g = z - y - h * f(tNext, z)
      const dg = 1 - h * partialY(f, tNext, z)
      // A near-zero denominator (h * ∂f/∂y ≈ 1) makes the correction g/dg
      // blow up — bail out and keep the current z rather than let a tiny,
      // noisy dg (from the numerical derivative) amplify into a huge,
      // wrong jump. 1e-8 is well above that noise floor.
      if (Math.abs(dg) < 1e-8) break
      const step = g / dg
      if (!Number.isFinite(step) || Math.abs(step) > 1e6 * Math.max(1, Math.abs(z))) break
      z = z - step
    }

    y = z
    t = tNext
    ts.push(t)
    ys.push(y)
  }
  return {ts, ys}
}

// 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, implicitEulerValue, trueValue,
//     eulerError, implicitEulerError }, 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}

  // Reference value: RK4 with a very fine step count, far finer than any n
  // tested below.
  const refSteps = 2000
  const trueValue = VM.numerical.rk4Solve(f, t0, y0, b, refSteps).ys[refSteps]

  const rows = []
  let n = 1
  while (n <= maxN) {
    const eulerValue = VM.numerical.eulerSolve(f, t0, y0, b, n).ys[n]
    const implicitEulerValue = implicitEulerSolve(f, t0, y0, b, n).ys[n]

    const eulerError = Math.abs(eulerValue - trueValue)
    const implicitEulerError = Math.abs(implicitEulerValue - trueValue)

    rows.push({n, eulerValue, implicitEulerValue, trueValue, eulerError, implicitEulerError})
    n++
  }

  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: [], implicitEulerTs: [], implicitEulerYs: [], 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 implicitEuler = implicitEulerSolve(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})

  // Cap at 1e8 so a diverging method (e.g. unstable explicit Euler) doesn't
  // blow up the axis range.
  const yValues = []
  for (const y of reference.ys) if (Number.isFinite(y) && Math.abs(y) <= 1e8) yValues.push(y)
  for (const y of euler.ys) if (Number.isFinite(y) && Math.abs(y) <= 1e8) yValues.push(y)
  for (const y of implicitEuler.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,
    implicitEulerTs: implicitEuler.ts, implicitEulerYs: implicitEuler.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, showEuler, showImplicitEuler) => {
    const traces = [
      // The fine reference solution curve
      {x: data.refTs, y: data.refYs, type: "scatter", mode: "lines", name: "Reference solution", showlegend: true, line: {color: "#2563eb", 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: true,
        line: {color: "#dc2626", width: 2}, marker: {color: "#dc2626", size: 6}
      })
    }
    if (showImplicitEuler && data.implicitEulerTs.length > 0) {
      traces.push({
        x: data.implicitEulerTs, y: data.implicitEulerYs, type: "scatter", mode: "lines+markers", name: "Implicit Euler", showlegend: true,
        line: {color: "#16a34a", width: 2}, marker: {color: "#16a34a", size: 6}
      })
    }

    const layout = {
      margin: {l: 0, r: 0, t: 0, b: 0},
      xaxis: {title: "t", range: [data.tRange.lo, data.tRange.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, showEuler, showImplicitEuler)

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 implicitEulerName = "Implicit Euler"
  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: "implicitEulerValue", stroke: () => implicitEulerName}),
    Plot.dot(result.rows, {x: "n", y: "implicitEulerValue", fill: () => implicitEulerName, r: 3})
  ]
  const colorDomain = [eulerName, implicitEulerName, trueName]
  const colorRange = ["#dc2626", "#16a34a", "#111827"]

  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 both 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 eulerPts = []
  const implicitEulerPts = []
  for (const row of result.rows) {
    if (Number.isFinite(row.eulerError) && row.eulerError > 0) {
      eulerPts.push({x: Math.log10(row.n), y: Math.log10(row.eulerError)})
    }
    if (Number.isFinite(row.implicitEulerError) && row.implicitEulerError > 0) {
      implicitEulerPts.push({x: Math.log10(row.n), y: Math.log10(row.implicitEulerError)})
    }
  }

  const eulerName = "Euler error"
  const implicitEulerName = "Implicit Euler error"
  const marks = [
    Plot.dot(eulerPts, {x: "x", y: "y", fill: () => eulerName, r: 4}),
    Plot.dot(implicitEulerPts, {x: "x", y: "y", fill: () => implicitEulerName, r: 4})
  ]
  const colorDomain = [eulerName, implicitEulerName]
  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 eulerFit = VM.numerical.linearRegression(eulerPts)
  if (eulerFit) {
    const lineName = `Euler order ≈ ${(-eulerFit.slope).toFixed(2)}`
    marks.push(Plot.line(
      [{x: eulerFit.xlo, y: eulerFit.slope * eulerFit.xlo + eulerFit.intercept}, {x: eulerFit.xhi, y: eulerFit.slope * eulerFit.xhi + eulerFit.intercept}],
      {x: "x", y: "y", stroke: () => lineName, strokeDasharray: "6,4"}
    ))
    colorDomain.push(lineName)
    colorRange.push("#f87171")
    legendDomain.push(lineName)
    legendRange.push("#f87171")
  }

  const implicitEulerFit = VM.numerical.linearRegression(implicitEulerPts)
  if (implicitEulerFit) {
    const lineName = `Implicit Euler order ≈ ${(-implicitEulerFit.slope).toFixed(2)}`
    marks.push(Plot.line(
      [{x: implicitEulerFit.xlo, y: implicitEulerFit.slope * implicitEulerFit.xlo + implicitEulerFit.intercept}, {x: implicitEulerFit.xhi, y: implicitEulerFit.slope * implicitEulerFit.xhi + implicitEulerFit.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
}

NoteValue table

// Value table — shows the entire run, not just the current n, with both
// 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.implicitEulerValue.toPrecision(10),
      row.trueValue.toPrecision(10),
      row.eulerError.toExponential(4),
      row.implicitEulerError.toExponential(4)
    ])
  }

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