Resources Numerical accuracy Calculation policy

Solvers shared by the process engine

The linear solve, Newton iteration, ODE integrator and bracketed root finder that the flowsheet, reactors and design specifications call.

How the precision is kept

All four work on plain arrays of 64-bit numbers. Their results are as close as their stopping rules: an absolute residual norm for Newton, relative and absolute error per step for the integrator, an absolute function value or a relative bracket for the root finder. Callers pass their own tolerances.

Open item. Return a flag from rk45 and illinois when the step or iteration limit is reached before the tolerance is met.

Calculations under this policy

4 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.

process/numeric.js CadCraft 0.8.2 64-bit 4 calculations Reviewed Nov 5, 2026

Tests tools/process-test/test.mjs

Linear solve

64-bit

x from A x = b for a small dense system.

Gaussian elimination with row exchange on the largest pivot, then back substitution

Partial pivoting on the absolute value in the column; the rows are not scaled first. Only a pivot below 1e-300 is called singular, so an ill-conditioned system is solved as it stands with no estimate of the digits lost.

process/numeric.js
lines 4 to 16
CadCraft 0.8.2

/** Solve A x = b (A n×n row arrays) by Gaussian elimination with partial pivoting; returns null if singular. */
export function solveLinear(A, b) {
  const n = b.length, M = A.map((r, i) => [...r, b[i]]);
  for (let c = 0; c < n; c++) {
    let p = c; for (let r = c + 1; r < n; r++) if (Math.abs(M[r][c]) > Math.abs(M[p][c])) p = r;
    if (Math.abs(M[p][c]) < 1e-300) return null;
    [M[c], M[p]] = [M[p], M[c]];
    for (let r = c + 1; r < n; r++) { const f = M[r][c] / M[c][c]; if (f) for (let k = c; k <= n; k++) M[r][k] -= f * M[c][k]; }
  }
  const x = new Array(n);
  for (let i = n - 1; i >= 0; i--) { let s = M[i][n]; for (let k = i + 1; k < n; k++) s -= M[i][k] * x[k]; x[i] = s / M[i][i]; }
  return x;
}

Newton iteration

64-bit

The solution of a set of nonlinear equations F(x) = 0.

J Δx = −F, Jᵢⱼ ≈ (Fᵢ(x + h eⱼ) − Fᵢ(x)) ÷ h, h = 1e-7 × max(|xⱼ|, 1e-6); x ← x + λ Δx, λ halved until |F| falls

The Jacobian is a one-sided difference, good to about half the digits, which slows convergence but does not limit the answer: the stop is on the residual norm, an absolute number. When no step lowers the residual the result is still reported as converged if the norm is below 1000 × tol.

process/numeric.js
lines 17 to 34
CadCraft 0.8.2

/** Newton on F(x)=0 with a forward-difference Jacobian, damping and an optional feasibility test. → { x, ok, iter, norm } */
export function newtonFD(F, x0, { tol = 1e-10, maxIter = 100, feasible = null, rel = 1e-7 } = {}) {
  let x = x0.slice(), f = F(x), nrm = Math.hypot(...f);
  for (let it = 0; it < maxIter; it++) {
    if (nrm < tol) return { x, ok: true, iter: it, norm: nrm };
    const n = x.length, J = Array.from({ length: f.length }, () => new Array(n));
    for (let j = 0; j < n; j++) { const h = rel * Math.max(Math.abs(x[j]), 1e-6), xp = x.slice(); xp[j] += h; const fp = F(xp); for (let i = 0; i < f.length; i++) J[i][j] = (fp[i] - f[i]) / h; }
    const dx = solveLinear(J, f.map((v) => -v)); if (!dx) return { x, ok: false, iter: it, norm: nrm };
    let lam = 1, moved = false;
    for (let k = 0; k < 40; k++) {
      const xn = x.map((v, i) => v + lam * dx[i]);
      if (!feasible || feasible(xn)) { const fn = F(xn), nn = Math.hypot(...fn); if (nn < nrm || nn < tol) { x = xn; f = fn; nrm = nn; moved = true; break; } }
      lam *= 0.5;
    }
    if (!moved) return { x, ok: nrm < tol * 1e3, iter: it, norm: nrm };
  }
  return { x, ok: nrm < tol * 1e3, iter: maxIter, norm: nrm };
}

Adaptive Runge–Kutta integration

64-bit

y(x₁) for a system of ordinary differential equations, with the accepted points along the way.

The local error per step is held to the tolerances; the global error is not estimated. A step shorter than 1e-14 × max(1, |x₁|) is accepted whatever its error. If 200000 steps are used up, the value at the last accepted point is returned with no flag.

process/numeric.js
lines 35 to 49
CadCraft 0.8.2

/** Dormand–Prince RK45, adaptive; f(x, y) → dy. Returns { y, steps, points[] } with the accepted points (x, y). */
export function rk45(f, y0, x0, x1, { rtol = 1e-9, atol = 1e-12, h0 = null, maxSteps = 200000 } = {}) {
  const c = [0, 1 / 5, 3 / 10, 4 / 5, 8 / 9, 1, 1];
  const a = [[], [1 / 5], [3 / 40, 9 / 40], [44 / 45, -56 / 15, 32 / 9], [19372 / 6561, -25360 / 2187, 64448 / 6561, -212 / 729], [9017 / 3168, -355 / 33, 46732 / 5247, 49 / 176, -5103 / 18656], [35 / 384, 0, 500 / 1113, 125 / 192, -2187 / 6784, 11 / 84]];
  const b5 = [35 / 384, 0, 500 / 1113, 125 / 192, -2187 / 6784, 11 / 84, 0], b4 = [5179 / 57600, 0, 7571 / 16695, 393 / 640, -92097 / 339200, 187 / 2100, 1 / 40];
  let x = x0, y = y0.slice(), h = h0 || (x1 - x0) / 100, steps = 0; const points = [[x, y.slice()]];
  while (x < x1 && steps < maxSteps) {
    if (x + h > x1) h = x1 - x;
    const k = []; for (let s = 0; s < 7; s++) { const ys = y.map((v, i) => { let t = 0; for (let j = 0; j < s; j++) t += a[s][j] * k[j][i]; return v + h * t; }); k.push(f(x + c[s] * h, ys)); }
    let err = 0; const yn = y.map((v, i) => { let t = 0, e = 0; for (let s = 0; s < 7; s++) { t += b5[s] * k[s][i]; e += (b5[s] - b4[s]) * k[s][i]; } const sc = atol + rtol * Math.max(Math.abs(v), Math.abs(v + h * t)); err = Math.max(err, Math.abs(h * e) / sc); return v + h * t; });
    if (err <= 1 || h < 1e-14 * Math.max(1, Math.abs(x1))) { x += h; y = yn; points.push([x, y.slice()]); }
    h *= Math.min(5, Math.max(0.2, 0.9 * Math.pow(err || 1e-10, -0.2))); steps++;
  }
  return { y, steps, points };
}

Bracketed root

64-bit

A root of one equation between two points where the function changes sign.

c = (f_a b − f_b a) ÷ (f_a − f_b); the retained end's value is halved when the same side is kept twice

Stops when |f(c)| is below tol (an absolute value, 1e-12 by default) or the bracket is narrower than tol × (1 + |c|). After 200 iterations the last point is returned with no flag. No sign change at the ends returns null.

process/numeric.js
lines 50 to 60
CadCraft 0.8.2

/** Bracketed root find on [a,b] (f(a), f(b) of opposite sign) by bisection+secant (Illinois); deterministic. */
export function illinois(f, a, b, { tol = 1e-12, maxIter = 200 } = {}) {
  let fa = f(a), fb = f(b), side = 0; if (fa === 0) return a; if (fb === 0) return b; if (fa * fb > 0) return null;
  let c = a;
  for (let i = 0; i < maxIter; i++) {
    c = (fa * b - fb * a) / (fa - fb); const fc = f(c);
    if (Math.abs(fc) < tol || Math.abs(b - a) < tol * (1 + Math.abs(c))) return c;
    if (fc * fb > 0) { b = c; fb = fc; if (side === -1) fa /= 2; side = -1; } else { a = c; fa = fc; if (side === 1) fb /= 2; side = 1; }
  }
  return c;
}

In detail

Method: Gaussian elimination with partial pivoting; Newton with a forward-difference Jacobian and step halving; Dormand–Prince 5(4) with step control; the Illinois variant of false position.
Number types: 64-bit numbers in plain arrays. The Dormand–Prince coefficients are written as fractions and evaluated in 64-bit.
Tolerances: solveLinear: a pivot below 1e-300 (absolute) means singular. newtonFD: residual norm below tol (1e-10 by default), Jacobian step 1e-7 × max(|x|, 1e-6), up to 40 halvings, 100 iterations. rk45: rtol 1e-9, atol 1e-12 by default, 200000 steps. illinois: |f| < 1e-12 or bracket below 1e-12 × (1 + |c|), 200 iterations.
Checked by: tools/process-test/test.mjs compares plug-flow and stirred-tank reactors, which call rk45 and newtonFD, with their closed-form conversions; it runs the library copy of this file, which is byte-identical.
Limits: rk45 and illinois return their last value without a flag when the step or iteration limit is reached. Newton's tolerance is on the unscaled residual, so the caller has to scale its equations.

Approach: 64-bit, stated limit

The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.

Every policy that takes this approach