Operating point and response over time of electrical, thermal, hydraulic and mechanical networks built from lumped elements.
How the precision is kept
Each time step is a Newton solve in 64-bit with a finite-difference Jacobian. How close the response is over time is set by the step: fixed at a thousandth of the run by default, or controlled by step doubling to a relative 1e-6. Newton's own stop includes a floor of 1e-11 of the largest unknown, so small unknowns beside large ones are found less closely.
Open item. Make the reported diode current follow the same straight-line law above 40 n Vt that the solved equations use.
Calculations under this policy
12 calculations in 3 source files. 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.
model.jsShared shell, every app64-bit4 calculationsReviewed
The flow through a resistor, damper, orifice, diode or user law for a given difference across it.
resistor: i = v ÷ R; orifice: i = v ÷ √(K √(v² + ε²)); diode: i = Is (e^(v ÷ nVt) − 1), a straight line above 40 nVt
The orifice law is smoothed with ε = 1e-3 (in the unit of the across variable) so its slope is finite at zero; well above ε it equals √(|v| ÷ K). The diode uses exp(x) − 1 and continues as a straight line above 40 n Vt. Dampers multiply by the coefficient; resistors divide.
model.js
lines 96 to 103
case 'R': case 'D': case 'N': case 'O': {
const v = V(x, p.a) - V(x, p.b); let i;
if (p.kind === 'R') { const par = val(p, p.T.p, t); i = p.T.conductance ? par * v : v / par; }
else if (p.kind === 'O') { const K = val(p, 'K', t), eps = e.eps ?? 1e-3; i = v / Math.sqrt(K * Math.sqrt(v * v + eps * eps)); }
else if (p.kind === 'D') { const nv = e.n * e.Vt, lim = 40 * nv; i = v > lim ? e.Is * (Math.exp(40) - 1) * (1 + (v - lim) / nv) : e.Is * (Math.exp(v / nv) - 1); }
else i = e.f(v, t);
if (p.a >= 0) r[p.a] += i; if (p.b >= 0) r[p.b] -= i; break;
}
Storage elements
64-bit
Capacitors, masses, tanks, inductors and springs: what they store and how it enters the equations.
capacitive: i + d(C v)/dt = 0; inductive: v + d(L i)/dt = 0; spring: L = 1 ÷ k
Stored quantities are C·v and L·i in 64-bit. A spring's compliance is the reciprocal of its stiffness, one extra rounding. At the DC operating point capacitive branches carry no flow and inductive ones no drop.
model.js
lines 106 to 121
case 'C': {
const i = x[p.br], v = V(x, p.a) - V(x, p.b); if (p.a >= 0) r[p.a] += i; if (p.b >= 0) r[p.b] -= i;
const Cv = typeof e.q === 'function' ? e.q(v) : val(p, p.T.p, t) * v;
if (mode === 'ic' && e.ic != null) r[p.br] += v - e.ic;
else if (mode !== 'tran') r[p.br] += i;
else { r[p.br] += i; S[p.br] -= Cv; }
break;
}
case 'L': {
const i = x[p.br], v = V(x, p.a) - V(x, p.b), Lp = val(p, p.T.p, t), L = p.T.inverse ? 1 / Lp : Lp;
if (p.a >= 0) r[p.a] += i; if (p.b >= 0) r[p.b] -= i;
if (mode === 'ic' && e.ic != null) r[p.br] += i - e.ic;
else if (mode !== 'tran') r[p.br] += v;
else { r[p.br] += v; S[p.br] -= L * i; }
break;
}
Reported flow through an element
64-bit
The current, heat flow, flow or force reported for an element.
the element's law at the solved state; diode: Is (e^(min(v, 40 nVt) ÷ nVt) − 1)
For a diode above 40 n Vt this function returns the constant Is (e⁴⁰ − 1), while the equations that were solved used the straight-line continuation, so the reported current is lower than the one in the balance. Below that voltage the two agree.
model.js
lines 135 to 147
/** through value of element `id` (from a to b; for a transformer the current into port 1) at state x, time t */
through(id, x, t) {
const p = parts.find((q) => q.e.id === id); if (!p) return NaN; const { e } = p;
const v = p.kind === 'X' || p.kind === 'Y' ? 0 : V(x, p.a) - V(x, p.b);
switch (p.kind) {
case 'R': { const par = val(p, p.T.p, t); return p.T.conductance ? par * v : v / par; }
case 'O': { const K = val(p, 'K', t), eps = e.eps ?? 1e-3; return v / Math.sqrt(K * Math.sqrt(v * v + eps * eps)); }
case 'D': { const nv = e.n * e.Vt; return e.Is * (Math.exp(Math.min(v, 40 * nv) / nv) - 1); }
case 'N': return e.f(v, t);
case 'J': return val(p, p.T.p, t);
default: return x[p.br];
}
},
solve.jsShared shell, every app64-bit6 calculationsReviewed
Singular when a pivot is at most 1e-13 of the largest matrix entry, so the test scales with the matrix. No refinement step.
solve.js
lines 20 to 33
function solveLinear(A, b) {
const n = b.length, M = A.map((r, i) => [...r, b[i]]);
let scale = 0; for (const r of A) for (const v of r) scale = Math.max(scale, Math.abs(v));
const tiny = Math.max(scale, 1e-300) * 1e-13;
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]) <= tiny) 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
Finds the state that satisfies the equations at one time.
J_ij = (F_i(x + h e_j) − F_i(x − h e_j)) ÷ 2h, h = 1e-7 × max(|x_j|, 1); stop when |Δxᵢ| ≤ rtol |xᵢ| + atol + 1e-11 × max(1, largest |x|)
The Jacobian is a central difference; it steers the iteration and does not move the root. The stop adds a floor of 1e-11 times the largest unknown to every component, so an unknown a million times smaller than the largest is only found to about 5 digits. A stalled line search is accepted if the step is already below that stop with a hundred times the floor, if the residual fell by a factor of 1e7, or if it is below 1e-12.
solve.js
lines 35 to 59
/** Newton with a central-difference Jacobian and backtracking. F(x) → residual array. → { x, ok, singular, iter } */
export function newton(F, x0, { rtol = 1e-9, atol = 1e-12, maxIter = 60 } = {}) {
let x = x0.slice(), f = F(x), nrm = Math.hypot(...f);
const n = x.length, nrm0 = nrm;
for (let it = 0; it < maxIter; it++) {
if (!Number.isFinite(nrm)) return { x, ok: false, iter: it };
const J = Array.from({ length: n }, () => new Array(n));
for (let j = 0; j < n; j++) {
const h = 1e-7 * Math.max(Math.abs(x[j]), 1), xp = x.slice(), xm = x.slice(); xp[j] += h; xm[j] -= h;
const fp = F(xp), fm = F(xm); for (let i = 0; i < n; i++) J[i][j] = (fp[i] - fm[i]) / (2 * h);
}
const dx = solveLinear(J, f.map((v) => -v)); if (!dx) return { x, ok: false, singular: true, iter: it };
let lam = 1, accepted = false, xn, fn, nn;
for (let k = 0; k < 30; k++) { xn = x.map((v, i) => v + lam * dx[i]); fn = F(xn); nn = Math.hypot(...fn); if (nn < nrm || nn < 1e-14) { accepted = true; break; } lam /= 2; }
const sc = 1e-11 * Math.max(1, ...x.map(Math.abs)); // roundoff floor of the whole vector: a component near zero cannot be asked for more
if (!accepted) { // stalled: accept if the step is already below the roundoff floor
const small = dx.every((d, i) => Math.abs(d) <= rtol * Math.abs(x[i]) + atol + sc * 100); return { x, ok: small || nrm <= 1e-7 * nrm0 || nrm < 1e-12, iter: it };
}
const conv = dx.every((d, i) => Math.abs(lam * d) <= rtol * Math.abs(xn[i]) + atol + sc);
x = xn; f = fn; nrm = nn;
if (conv) return { x, ok: true, iter: it + 1 };
if (nrm < 1e-15) return { x, ok: true, iter: it + 1 };
}
return { x, ok: false, iter: maxIter };
}
Operating point
64-bit
The steady state, or the state set by the given initial values.
one Newton solve from zero
Capacitive branches are open and inductive ones shorted, or the given initial values are imposed. Tolerances are 1e-9 relative and 1e-12 absolute plus the floor of the Newton stop.
solve.js
lines 83 to 90
export function operatingPoint(model, { init = 'dc', t0 = 0 } = {}) {
const c = model.compile(); if (!c.ok) return { ok: false, problems: c.errors.map((e, i) => problem('build.' + i, 'error', e, [])) };
const sys = c.system, n = sys.size;
const F = (x) => { const r = new Array(n).fill(0), S = new Array(n).fill(0); sys.eval(x, t0, init === 'ic' ? 'ic' : 'dc', r, S); return r; };
const s = newton(F, new Array(n).fill(0));
if (!s.ok) return { ok: false, problems: s.singular ? diagnose(sys, ' at the operating point') : [problem('noconverge', 'error', 'The operating point did not converge.', [{ kind: 'model', id: 'all' }])], system: sys };
return { ok: true, x: s.x, system: sys, problems: [], v: (name) => sys.across(name, s.x), i: (id) => sys.through(id, s.x, t0) };
}
One Newton solve per step with rtol 1e-10. The trapezoid carries the derivative by the recursion Ṡ = (2 ÷ h)(S − S₀) − Ṡ₀, in which an error in Ṡ₀ changes sign each step and does not die out. The first step after the start or a breakpoint is a short backward-Euler step.
solve.js
lines 103 to 110
const advance = (s, h, beOnce = false) => {
const useTrap = trap && !beOnce, t1 = s.t + h, c = useTrap ? 2 / h : 1 / h;
const F = (x) => { const { r, S } = parts(x, t1); for (const i of dynRows) r[i] += c * (S[i] - s.S[i]) - (useTrap ? s.Sd[i] : 0); return r; };
const sol = newton(F, s.x, { rtol: 1e-10, atol: 1e-13 });
if (!sol.ok) return sol.singular ? { singular: true } : null;
const { S } = parts(sol.x, t1), Sd = zero(); for (const i of dynRows) Sd[i] = c * (S[i] - s.S[i]) - (useTrap ? s.Sd[i] : 0);
return { t: t1, x: sol.x, S, Sd };
};
Step-size control
64-bit
Chooses the time step from an estimate of the error, when adaptive stepping is on.
err = max |x(two half steps) − x(one step)| ÷ (rtol |x| + atol), ÷ 3 for the trapezoid; h ← h × clamp(0.9 err^(−1/3), 0.2, 5)
The estimate compares one full step with two half steps; the two-half-step result is kept without extrapolation. Defaults are rtol 1e-6 and atol 1e-9. A step at the minimum size is accepted whatever its error. Without adaptive stepping the step is fixed and no error is estimated.
solve.js
lines 127 to 137
// step doubling
const a = advance(st, hh); const b1 = a && !a.singular ? advance(st, hh / 2) : null; const b2 = b1 && !b1.singular ? advance(b1, hh / 2) : null;
if (!a || !b2 || a.singular) {
if (hh <= hMin * 2) { problems.push(problem('noconverge', 'error', `Newton did not converge at t = ${(st.t + hh).toPrecision(5)} s even at the smallest step.`, [{ kind: 'model', id: 'all' }])); break; }
h = hh / 4; rejected++; continue;
}
let err = 0; for (let i = 0; i < n; i++) err = Math.max(err, Math.abs(b2.x[i] - a.x[i]) / (rtol * Math.abs(b2.x[i]) + atol));
err /= trap ? 3 : 1;
if (err <= 1 || hh <= hMin * 2) { st = b2; steps++; T.push(st.t); X.push(st.x); if (hitBp && st.t >= nb - 1e-15) restart = true; } else rejected++;
const fac = err > 0 ? Math.min(5, Math.max(0.2, 0.9 * Math.pow(err, trap ? -1 / 3 : -1 / 2))) : 5;
h = Math.min(hMax, Math.max(hMin, hh * fac));
Value between stored times
64-bit
A node value at a time between two computed points.
y = y₀ + (y₁ − y₀)(t − t₀) ÷ (t₁ − t₀)
Values between stored points are straight lines, which is a lower order than the trapezoidal rule that produced the points.
solve.js
lines 151 to 155
function interp(T, y, t) {
if (t <= T[0]) return y[0]; if (t >= T[T.length - 1]) return y[y.length - 1];
let lo = 0, hi = T.length - 1; while (hi - lo > 1) { const m = (lo + hi) >> 1; (T[m] <= t ? (lo = m) : (hi = m)); }
return y[lo] + (y[hi] - y[lo]) * (t - T[lo]) / (T[hi] - T[lo]);
}
modelica.jsShared shell, every app64-bit2 calculationsReviewed
The value of a parameter written as an expression, such as 2*3.3e3.
numbers, parameters, + − × ÷ ^, brackets, pi
parseFloat for numbers and 64-bit operators; ^ is Math.pow and groups to the right. A leading minus applies to the whole power that follows, so −2^2 is −4.
modelica.js
lines 18 to 35
function evalExpr(text, vars) {
const toks = text.match(/\d+\.?\d*(?:[eE][+-]?\d+)?|\.\d+(?:[eE][+-]?\d+)?|[A-Za-z_]\w*|[()+\-*/^]/g) || []; let i = 0;
const peek = () => toks[i], next = () => toks[i++];
function prim() {
const t = next();
if (t === undefined) throw new Error('unexpected end of expression');
if (t === '(') { const v = add(); if (next() !== ')') throw new Error('missing )'); return v; }
if (t === '-') return -pow(); if (t === '+') return pow();
if (/^[\d.]/.test(t)) return parseFloat(t);
if (Object.prototype.hasOwnProperty.call(vars, t)) return vars[t];
if (t === 'pi') return Math.PI;
throw new Error('unknown name ' + t);
}
function pow() { const b = prim(); if (peek() === '^') { next(); return Math.pow(b, pow()); } return b; }
function mul() { let v = pow(); while (peek() === '*' || peek() === '/') v = next() === '*' ? v * pow() : v / pow(); return v; }
function add() { let v = mul(); while (peek() === '+' || peek() === '-') v = next() === '+' ? v + mul() : v - mul(); return v; }
const v = add(); if (i < toks.length) throw new Error('unexpected ' + toks[i]); return v;
}
Thermal conductance to resistance
64-bit
A thermal conductor from a model file as the engine's thermal resistance.
R = 1 ÷ G
The conductance is turned into a resistance by one division; the model then divides the temperature difference by it, two roundings where a multiplication by G would be one. A fixed temperature defaults to 293.15 K.
Differential-algebraic form g(x, t) + dS/dt = 0; trapezoidal rule or backward Euler; Newton's method with a central-difference Jacobian and backtracking; step-doubling (Richardson) error estimate; a reader for a subset of the Modelica language.
Number types
64-bit numbers in plain arrays.
Tolerances
Newton: |Δx| ≤ rtol |x| + atol + 1e-11 × max(1, largest |x|), with rtol 1e-9 and atol 1e-12 at the operating point and 1e-10 and 1e-13 in a time step; 60 iterations. Jacobian step 1e-7 × max(|x|, 1). Linear solve singular at a pivot of at most 1e-13 of the largest entry. Adaptive stepping: rtol 1e-6, atol 1e-9. Orifice smoothing 1e-3. Diode law linear beyond 40 n Vt.
Checked by
The engines test compares a cooling example at 100 s with 293.15 + 80 e⁻¹ within 1e-3, and a divider midpoint with 6 V within 1e-6.
Limits
The diode's default thermal voltage is 0.025852 V. The current reported for a diode above 40 n Vt is held at Is (e⁴⁰ − 1), which is not the current the solver balanced there. Values between stored times are linear interpolation.
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.