Butcher Tableau Explorer

// VM is the shared utility library, loaded globally via _includes/head-scripts.html (js/**). This page uses VM.expressions.makeNumber, VM.expressions.makeRational.
VM = window.VM

viewof stagesInput = {
  const params = new URLSearchParams(window.location.search)
  const input = Inputs.text({label: "Number of stages (s)", value: params.get("s") ?? "2", placeholder: "1-8"})
  input.classList.add("ojs-auto")
  return input
}

// stages clamps the stage-count field to a usable integer in [1, 8] — large
// enough for common methods (RK4 = 4, Fehlberg = 6, Dormand-Prince = 7)
// without letting the tableau or the Taylor-series computation below grow
// unreasonably large.
stages = {
  const math = window.math
  const parsed = VM.expressions.makeNumber(math, stagesInput)
  let n = 2
  if (parsed !== null && Number.isFinite(parsed)) n = Math.round(parsed)
  if (n < 1) n = 1
  if (n > 8) n = 8
  return n
}

// rat: minimal exact-rational arithmetic used by the symbolic engine below.
// Every {num, den} value is kept reduced with den > 0.
rat = {
  const gcd = (a, b) => {
    a = Math.abs(a)
    b = Math.abs(b)
    while (b !== 0) {
      const remainder = a % b
      a = b
      b = remainder
    }
    if (a === 0) return 1
    return a
  }

  const reduce = (num, den) => {
    if (den < 0) {
      num = -num
      den = -den
    }
    const divisor = gcd(num, den)
    return {num: num / divisor, den: den / divisor}
  }

  const of = n => ({num: n, den: 1})
  const add = (a, b) => reduce(a.num * b.den + b.num * a.den, a.den * b.den)
  const mul = (a, b) => reduce(a.num * b.num, a.den * b.den)
  const div = (a, b) => reduce(a.num * b.den, a.den * b.num)
  const isZero = a => a.num === 0

  const toPlain = a => {
    if (a.den === 1) return `${a.num}`
    return `${a.num}/${a.den}`
  }

  const toLatex = a => {
    let sign = ""
    let num = a.num
    if (num < 0) {
      sign = "-"
      num = -num
    }
    if (a.den === 1) return `${sign}${num}`
    return `${sign}\\frac{${num}}{${a.den}}`
  }

  return {of, add, mul, div, isZero, toPlain, toLatex}
}

// formatLinearCombination(coeffs, symbolNames) renders coeffs[0]*symbolNames[0]
// + coeffs[1]*symbolNames[1] + ... as LaTeX, dropping zero terms and writing
// coefficients of exactly 1 as just the symbol — or "" if every coefficient
// is zero.
formatLinearCombination = (coeffs, symbolNames) => {
  let str = ""
  for (let idx = 0; idx < coeffs.length; idx++) {
    const coeff = coeffs[idx]
    if (rat.isZero(coeff)) continue

    let sign = "+"
    let magnitude = coeff
    if (coeff.num < 0) {
      sign = "-"
      magnitude = {num: -coeff.num, den: coeff.den}
    }
    const magLatex = rat.toLatex(magnitude)
    let piece = `${magLatex}${symbolNames[idx]}`
    if (magLatex === "1") piece = symbolNames[idx]

    if (str === "") {
      if (sign === "-") str = `-${piece}`
      else str = piece
    } else {
      str += ` ${sign} ${piece}`
    }
  }
  return str
}

// viewof tableauInput builds an editable s x s matrix of a_ij inputs, plus a
// column of c_i inputs and a row of b_i inputs, rebuilding from all-zero
// entries whenever the stage count changes. A column to the right of each
// row shows the live k_i (or y_{i+1}) equation that row's current entries
// construct — c_i doesn't feed into it, since it never enters the
// autonomous computation.
viewof tableauInput = {
  const s = stages

  // Seed from the "a" (rows joined by ";", entries within a row joined by
  // ",") / "b" / "c" URL params when present, falling back to "0" for any
  // entry the URL doesn't cover (missing param, or a saved shape smaller
  // than the current stage count).
  const params = new URLSearchParams(window.location.search)
  const aParam = params.get("a")
  const bParam = params.get("b")
  const cParam = params.get("c")
  const aRows = aParam === null ? [] : aParam.split(";")
  const bVals = bParam === null ? [] : bParam.split(",")
  const cVals = cParam === null ? [] : cParam.split(",")

  const state = {a: [], b: [], c: []}
  for (let i = 0; i < s; i++) {
    const row = []
    const savedRow = i < aRows.length ? aRows[i].split(",") : []
    for (let j = 0; j < s; j++) {
      if (j < savedRow.length) row.push(savedRow[j])
      else row.push("0")
    }
    state.a.push(row)
  }
  for (let j = 0; j < s; j++) {
    if (j < bVals.length) state.b.push(bVals[j])
    else state.b.push("0")
  }
  for (let i = 0; i < s; i++) {
    if (i < cVals.length) state.c.push(cVals[i])
    else state.c.push("0")
  }

  const container = document.createElement("div")
  container.className = "butcher-tableau-wrap"

  const table = document.createElement("table")
  table.className = "butcher-tableau"
  const tbody = document.createElement("tbody")

  const emit = () => {
    container.value = {a: state.a, b: state.b, c: state.c}
    container.dispatchEvent(new CustomEvent("input"))
  }

  const makeCoeffInput = (initialValue, onChange) => {
    const input = document.createElement("input")
    input.type = "text"
    input.value = initialValue
    input.addEventListener("input", () => {
      const parsed = VM.expressions.makeRational(input.value)
      if (parsed === null) input.classList.add("bt-invalid")
      else input.classList.remove("bt-invalid")
      onChange(input.value)
    })
    return input
  }

  const setEquationCell = (td, latex) => {
    while (td.firstChild) td.removeChild(td.firstChild)
    td.appendChild(tex`${latex}`)
  }

  const stageRowLatex = i => {
    const coeffs = []
    const symbolNames = []
    for (let j = 0; j < s; j++) {
      const parsed = VM.expressions.makeRational(state.a[i][j])
      if (parsed === null) coeffs.push(rat.of(0))
      else coeffs.push(parsed)
      symbolNames.push(`k_{${j + 1}}`)
    }
    const combo = formatLinearCombination(coeffs, symbolNames)
    if (combo === "") return `k_{${i + 1}} = f(y_i)`
    return `k_{${i + 1}} = f\\!\\left(y_i + h\\left(${combo}\\right)\\right)`
  }

  const bRowLatex = () => {
    const coeffs = []
    const symbolNames = []
    for (let j = 0; j < s; j++) {
      const parsed = VM.expressions.makeRational(state.b[j])
      if (parsed === null) coeffs.push(rat.of(0))
      else coeffs.push(parsed)
      symbolNames.push(`k_{${j + 1}}`)
    }
    const combo = formatLinearCombination(coeffs, symbolNames)
    if (combo === "") return "y_{i+1} = y_i"
    return `y_{i+1} = y_i + h\\left(${combo}\\right)`
  }

  for (let i = 0; i < s; i++) {
    const tr = document.createElement("tr")

    const cTd = document.createElement("td")
    cTd.className = "bt-c"
    const row = i
    const cInput = makeCoeffInput(state.c[row], value => {
      state.c[row] = value
      emit()
    })
    cTd.appendChild(cInput)
    tr.appendChild(cTd)

    for (let j = 0; j < s; j++) {
      const td = document.createElement("td")
      td.className = "bt-a"
      const col = j
      const input = makeCoeffInput(state.a[row][col], value => {
        state.a[row][col] = value
        setEquationCell(eqTd, stageRowLatex(row))
        emit()
      })
      td.appendChild(input)
      tr.appendChild(td)
    }

    const eqTd = document.createElement("td")
    eqTd.className = "bt-eq"
    setEquationCell(eqTd, stageRowLatex(i))
    tr.appendChild(eqTd)

    tbody.appendChild(tr)
  }

  const bTr = document.createElement("tr")
  bTr.className = "bt-brow"
  const blankTd = document.createElement("td")
  blankTd.className = "bt-c"
  bTr.appendChild(blankTd)
  for (let j = 0; j < s; j++) {
    const td = document.createElement("td")
    td.className = "bt-b"
    const col = j
    const input = makeCoeffInput(state.b[col], value => {
      state.b[col] = value
      setEquationCell(bEqTd, bRowLatex())
      emit()
    })
    td.appendChild(input)
    bTr.appendChild(td)
  }
  const bEqTd = document.createElement("td")
  bEqTd.className = "bt-eq"
  setEquationCell(bEqTd, bRowLatex())
  bTr.appendChild(bEqTd)
  tbody.appendChild(bTr)

  table.appendChild(tbody)
  container.appendChild(table)
  container.value = {a: state.a, b: state.b, c: state.c}
  return container
}

// tableau parses every cell's raw text into a Rational, falling back to 0
// for any cell that doesn't parse — the tableau always computes something
// rather than blocking on a mid-edit or invalid entry.
tableau = {
  const s = stages
  const a = []
  for (let i = 0; i < s; i++) {
    const row = []
    for (let j = 0; j < s; j++) {
      const parsed = VM.expressions.makeRational(tableauInput.a[i][j])
      if (parsed === null) row.push(rat.of(0))
      else row.push(parsed)
    }
    a.push(row)
  }
  const b = []
  for (let j = 0; j < s; j++) {
    const parsed = VM.expressions.makeRational(tableauInput.b[j])
    if (parsed === null) b.push(rat.of(0))
    else b.push(parsed)
  }

  // Mirror the raw (unparsed) tableau into the URL as "s" plus "a"/"b"/"c"
  // lists — "a" is rows joined by ";" with entries within a row joined by
  // ",", "b"/"c" are just "," joined — so a tableau can be shared/reloaded
  // via link, same as every other page's URL-synced inputs.
  const aRows = []
  for (let i = 0; i < s; i++) aRows.push(tableauInput.a[i].join(","))

  const params = new URLSearchParams(window.location.search)
  params.set("s", String(s))
  params.set("a", aRows.join(";"))
  params.set("b", tableauInput.b.join(","))
  params.set("c", tableauInput.c.join(","))
  history.replaceState(null, "", `${window.location.pathname}?${params}`)

  return {a, b}
}

// maxOrder caps how many powers of h the Taylor expansion below is computed
// to. It only needs to reach far enough to find two nonzero terms, which
// for almost any tableau happens by h^2 — 5 leaves generous headroom for
// tableaus deliberately built to cancel the first few terms.
maxOrder = 5

// polyOps: a "poly" is a list of {coeff, counts} terms, one per distinct
// monomial in the symbols f, f', f'', ... — counts[m] is the exponent of
// the m-th derivative of f (evaluated at y_i) in that monomial. Terms
// sharing a monomial are always merged, so a poly has at most one term per
// monomial.
polyOps = {
  const key = counts => counts.join(",")

  const zero = () => []

  const symbol = (m, size) => {
    const counts = new Array(size).fill(0)
    counts[m] = 1
    return [{coeff: rat.of(1), counts}]
  }

  const add = (p1, p2) => {
    const map = new Map()
    for (const term of p1) map.set(key(term.counts), {coeff: term.coeff, counts: term.counts})
    for (const term of p2) {
      const k = key(term.counts)
      const existing = map.get(k)
      if (existing !== undefined) existing.coeff = rat.add(existing.coeff, term.coeff)
      else map.set(k, {coeff: term.coeff, counts: term.counts})
    }
    const result = []
    for (const term of map.values()) {
      if (!rat.isZero(term.coeff)) result.push(term)
    }
    return result
  }

  const scale = (p, scalar) => {
    const result = []
    for (const term of p) {
      const coeff = rat.mul(term.coeff, scalar)
      if (!rat.isZero(coeff)) result.push({coeff, counts: term.counts})
    }
    return result
  }

  const mul = (p1, p2) => {
    let result = zero()
    for (const t1 of p1) {
      for (const t2 of p2) {
        const counts = new Array(t1.counts.length)
        for (let m = 0; m < counts.length; m++) counts[m] = t1.counts[m] + t2.counts[m]
        const coeff = rat.mul(t1.coeff, t2.coeff)
        result = add(result, [{coeff, counts}])
      }
    }
    return result
  }

  const isZero = p => p.length === 0

  return {zero, symbol, add, scale, mul, isZero}
}

// seriesOps: a "series" is an array of length maxOrder+1 of polys, where
// series[n] is the coefficient of h^n.
seriesOps = {
  const zero = () => {
    const result = []
    for (let n = 0; n <= maxOrder; n++) result.push(polyOps.zero())
    return result
  }

  const add = (s1, s2) => {
    const result = []
    for (let n = 0; n <= maxOrder; n++) result.push(polyOps.add(s1[n], s2[n]))
    return result
  }

  const scale = (s, scalar) => {
    const result = []
    for (let n = 0; n <= maxOrder; n++) result.push(polyOps.scale(s[n], scalar))
    return result
  }

  // Truncated Cauchy product: (s1 * s2)[n] = sum_{i+j=n} s1[i] * s2[j]
  const mul = (s1, s2) => {
    const result = []
    for (let n = 0; n <= maxOrder; n++) result.push(polyOps.zero())
    for (let i = 0; i <= maxOrder; i++) {
      for (let j = 0; j <= maxOrder - i; j++) {
        const term = polyOps.mul(s1[i], s2[j])
        result[i + j] = polyOps.add(result[i + j], term)
      }
    }
    return result
  }

  const one = () => {
    const result = zero()
    result[0] = [{coeff: rat.of(1), counts: new Array(maxOrder + 1).fill(0)}]
    return result
  }

  const power = (s, m) => {
    let result = one()
    for (let i = 0; i < m; i++) result = mul(result, s)
    return result
  }

  return {zero, add, scale, mul, one, power}
}

// kSeries[i] is the h-power series for stage value k_i(h), computed order
// by order from k_i(h) = f(y_i + h * sum_j a_ij k_j(h)):
//   k_i(h) = sum_{m>=0} (1/m!) f^(m)(y_i) * (h * S_i(h))^m,  S_i(h) = sum_j a_ij k_j(h)
// The coefficient of h^n in k_i only depends on the coefficients of
// h^0..h^(n-1) in every k_j (each term above contributes at order >= m>=1
// beyond order 0), so every order is solvable directly in sequence — even
// for an implicit tableau (a_ii != 0) — without solving any equation.
kSeries = {
  const s = stages
  const a = tableau.a

  const fSym = []
  for (let m = 0; m <= maxOrder; m++) fSym.push(polyOps.symbol(m, maxOrder + 1))

  const k = []
  for (let i = 0; i < s; i++) k.push(seriesOps.zero())

  // Order 0: k_i(0) = f(y_i) for every stage.
  for (let i = 0; i < s; i++) k[i][0] = fSym[0]

  for (let n = 1; n <= maxOrder; n++) {
    // S[i] = sum_j a_ij k_j(h), using only the already-known orders 0..n-1
    // of every k_j — orders n..maxOrder aren't filled in yet at this point,
    // so they're correctly zero placeholders rather than unknowns.
    const S = []
    for (let i = 0; i < s; i++) {
      let Si = seriesOps.zero()
      for (let j = 0; j < s; j++) Si = seriesOps.add(Si, seriesOps.scale(k[j], a[i][j]))
      S.push(Si)
    }

    let factorial = 1
    for (let m = 1; m <= n; m++) {
      factorial *= m
      const invFactorial = rat.div(rat.of(1), rat.of(factorial))
      for (let i = 0; i < s; i++) {
        const Sim = seriesOps.power(S[i], m)
        const term = polyOps.scale(polyOps.mul(fSym[m], Sim[n - m]), invFactorial)
        k[i][n] = polyOps.add(k[i][n], term)
      }
    }
  }

  return k
}

// updateSeries[n] is the coefficient of h^(n+1) in y_{i+1} - y_i = h * sum_i b_i k_i(h).
updateSeries = {
  const s = stages
  const b = tableau.b
  let C = seriesOps.zero()
  for (let i = 0; i < s; i++) C = seriesOps.add(C, seriesOps.scale(kSeries[i], b[i]))
  return C
}

// differentiatePoly(poly) returns d/dt of the elementary-differential
// polynomial poly along the true solution y' = f(y), via the chain rule:
// for each symbol f^(m) appearing in a term, decrementing its exponent by
// one and multiplying by f^(m+1) f contributes one summand — since
// d/dt f^(m)(y(t)) = f^(m+1)(y(t)) y'(t) = f^(m+1) f.
differentiatePoly = poly => {
  let result = polyOps.zero()
  for (const term of poly) {
    for (let m = 0; m < term.counts.length - 1; m++) {
      const exponent = term.counts[m]
      if (exponent === 0) continue
      const counts = term.counts.slice()
      counts[m] -= 1
      counts[m + 1] += 1
      counts[0] += 1
      const coeff = rat.mul(term.coeff, rat.of(exponent))
      result = polyOps.add(result, [{coeff, counts}])
    }
  }
  return result
}

// exactDerivatives[n] is the polynomial for y^(n+1)(t_i), the true
// solution's (n+1)-th derivative at the current point, for n = 0..maxOrder.
// exactDerivatives[0] = f (since y' = f); each later entry is the time
// derivative of the one before it.
exactDerivatives = {
  const derivatives = [polyOps.symbol(0, maxOrder + 1)]
  for (let n = 1; n <= maxOrder; n++) derivatives.push(differentiatePoly(derivatives[n - 1]))
  return derivatives
}

// exactSeries[n] is the coefficient of h^(n+1) in the true solution's own
// Taylor expansion y(t_i + h) - y_i = sum_n h^(n+1)/(n+1)! y^(n+1)(t_i).
exactSeries = {
  const result = []
  let factorial = 1
  for (let n = 0; n <= maxOrder; n++) {
    factorial *= (n + 1)
    const invFactorial = rat.div(rat.of(1), rat.of(factorial))
    result.push(polyOps.scale(exactDerivatives[n], invFactorial))
  }
  return result
}

// errorSeries[n] is the coefficient of h^(n+1) in y_{i+1} - y(t_i + h), the
// local truncation error of one step — the y_i term cancels between the two
// expansions above, so this is exactly updateSeries minus exactSeries.
errorSeries = {
  const result = []
  for (let n = 0; n <= maxOrder; n++) {
    result.push(polyOps.add(updateSeries[n], polyOps.scale(exactSeries[n], rat.of(-1))))
  }
  return result
}

// errorTerms: the leading (lowest) power of h (within [1, maxOrder]) whose
// coefficient in the local truncation error is not identically zero.
errorTerms = {
  const found = []
  for (let n = 0; n <= maxOrder; n++) {
    if (!polyOps.isZero(errorSeries[n])) found.push({power: n + 1, poly: errorSeries[n]})
    if (found.length === 1) break
  }
  return found
}

formatMonomial = counts => {
  const derivName = m => {
    if (m === 0) return "f"
    if (m === 1) return "f'"
    if (m === 2) return "f''"
    if (m === 3) return "f'''"
    return `f^{(${m})}`
  }

  let str = ""
  for (let m = 0; m < counts.length; m++) {
    const exponent = counts[m]
    if (exponent === 0) continue
    let factor = derivName(m)
    if (exponent > 1) factor = `${factor}^{${exponent}}`
    if (str !== "") str += " "
    str += factor
  }
  if (str === "") str = "1"
  return str
}

formatPoly = poly => {
  let str = ""
  for (const term of poly) {
    let sign = "+"
    let magnitude = term.coeff
    if (term.coeff.num < 0) {
      sign = "-"
      magnitude = {num: -term.coeff.num, den: term.coeff.den}
    }
    const monomialStr = formatMonomial(term.counts)
    const magLatex = rat.toLatex(magnitude)

    let piece = `${magLatex}\\,${monomialStr}`
    if (magLatex === "1" && monomialStr !== "1") piece = monomialStr

    if (str === "") {
      if (sign === "-") str = `-${piece}`
      else str = piece
    } else {
      str += ` ${sign} ${piece}`
    }
  }
  if (str === "") str = "0"
  return str
}

// errorLatex assembles "y_{i+1} - y(t_i+h) approx h^p(...) + O(h^n)" from
// the leading nonzero term of the local truncation error found above.
errorLatex = {
  if (errorTerms.length === 0) return `y_{i+1} - y(t_i + h) \\approx O(h^{${maxOrder + 1}})`

  const term = errorTerms[0]
  const polyLatex = formatPoly(term.poly)
  const remainderPower = term.power + 1

  return `y_{i+1} - y(t_i + h) \\approx h^{${term.power}}\\!\\left(${polyLatex}\\right) + O(h^{${remainderPower}})`
}

tex.block`${errorLatex}`

An \(s\)-stage Runge-Kutta method is specified by a Butcher tableau of nodes \(c_i\), coefficients \(a_{ij}\), and weights \(b_i\). Enter coefficients as integers, fractions, or decimals — each row’s panel shows the exact \(k_i\) (or \(y_{i+1}\)) equation your entries construct, and the box above computes the leading term of the local truncation error \(y_{i+1} - y(t_i+h)\) as a power series in \(h\), even for implicit tableaus. If that leading term is \(O(h^{p+1})\), the method is \(p\)-th order accurate. Try Euler’s method (\(s=1\), \(b_1=1\)), Heun’s method (\(s=2\), \(a_{21}=1\), \(b_1=b_2=1/2\)), or classical RK4 (\(s=4\), \(a_{21}=a_{32}=1/2\), \(a_{43}=1\), \(b_1=b_4=1/6\), \(b_2=b_3=1/3\)).