Module differentiate

Expand description

Numerical differentiation of a supplied function — finite differences, batched derivatives, and batched Jacobians, dispatched across ComputeBackend.

§This is NOT the FV gradient operator

crate::fv_operators implements the spatial finite-volume grad, div and laplacian over a mesh: they differentiate a field with respect to position, using face fluxes and cell volumes, and they are the right tool for a PDE discretisation. This module differentiates an arbitrary caller-supplied function with respect to its own arguments, by sampling it. If you are discretising a transport equation you want fv_operators; if you need df/dx of a closure, a property correlation or an ODE right-hand side, you are in the right place.

§The problem this exists to solve

crate::ode::OdeSystem::jacobian has a default body that is unimplemented!(). Any system that does not hand-code its Jacobian panics the moment crate::ode::Rosenbrock23 — the crate’s only stiff solver — asks for one, and inside crate::ode::parallel::integrate_ensemble that panic propagates out through the rayon scope. NumericalJacobian closes that hole: wrap a system, and Rosenbrock23 integrates it with no hand-coded Jacobian at all. Measured cost of doing so, on Van der Pol (mu = 5, x in [0, 10], tolerances 1e-8) — see “Measured cost against a hand-coded Jacobian” below — is 1.9x wall clock for a forward difference and 2.0x for a central one, for the same answer to all eight printed decimals.

§Provenance — a generalisation of two settled workspace conventions

Nothing here is a new algorithm. Both halves of the formulation are taken from implementations already working in this workspace, and the divergences are stated rather than left to be discovered.

The Jacobian assembly generalises:

crates/outram-park-fork-dwsim-libs/src/columns/linalg.rs:183
    pub fn finite_difference_jacobian<F>(f: &mut F, x: &[f64], epsilon: f64)
        -> Option<Array2<f64>>

itself a port of DWSIM’s NewtonRaphson.vb:669-705 (FunctionGradient), used there for the Naphtali-Sandholm column solver’s initial Broyden Jacobian. Kept from it: the central stencil, the relative perturbation, and the failure-is-Option convention — a caller never receives a matrix it cannot tell apart from a good one. dwsim-libs cannot be depended on from here (this crate has no internal workspace dependencies, by policy), so this is a reuse of formulation, not of code.

Three deliberate divergences from it, each because the alternative is a known defect:

This module	`finite_difference_jacobian`	Why
`h = rel * max(\|x\|, min_scale)`	`x(1±eps)`, or `eps` and `2eps` when `x == 0`	The `x == 0` branch silently switches to a one-sided stencil at a different step, so the scheme changes with the data. The `max` floor keeps one scheme everywhere.
A failed entry is `NaN` and the status says so	a non-finite entry is written as `0.0`	A zero is a plausible Jacobian entry. It cannot be detected downstream, and a Newton or Rosenbrock step built on it returns a wrong answer instead of an error.
Divides by the realised step `xp - xm`	divides by the requested `2epsx`	`x + h` is not representable, so the requested step is not the step taken. See `derivative`.

The step-size rule is the one already in:

crates/outram-park-fork-offbeat/src/rheology/aster/integration.rs:298
    pub fn newton_perturbed(...)   // h = perturbation * x.abs().max(1.0)
    pub fn perturbed_default() -> f64   // f64::EPSILON.cbrt()

which is upstream Code_Aster’s NEWTON_PERT. This module adopts both the max(|x|, 1) floor and eps^(1/3) for the central scheme verbatim — see CBRT_EPSILON and DiffSettings::step_for — and extends the same reasoning to the other three orders.

A third, narrower precedent — tampines-steam-tables’ w_ps_eqm_region4_finite_diff_vol (region_4_vap_liq_equilibrium/speed_of_sound_eqm.rs:83) — takes dv/dp|_s by central differences with a hard-coded dp = 1e-4 * p, clamped at the minimum table pressure. That is a relative step of 1e-4, sixteen times coarser than CBRT_EPSILON, which is the right call there because the IF97 flash it differentiates is far noisier than machine epsilon. It is recorded here as the standing reminder that the optimal step assumes the function is evaluated to rounding accuracy, and a caller whose function is noisier should raise DiffSettings::relative_step accordingly.

§Achievable accuracy is `sqrt(eps)` to `eps^(4/5)`, NEVER `eps`

This is the single most common misunderstanding about finite differences, so it is measured rather than asserted. Truncation error falls as h^p while round-off grows as eps/h; their sum is minimised at h ~ eps^(1/(p+1)), where the achievable accuracy is ~ eps^(p/(p+1)):

Scheme	Order `p`	Optimal `h`	Predicted accuracy
`DiffScheme::Forward` / `DiffScheme::Backward`	1	`sqrt(eps) = 1.49e-8`	`sqrt(eps) = 1.49e-8`
`DiffScheme::Central`	2	`eps^(1/3) = 6.06e-6`	`eps^(2/3) = 3.67e-11`
`DiffScheme::Central4th`	4	`eps^(1/5) = 7.40e-4`	`eps^(4/5) = 3.00e-13`

A caller expecting 1e-15 from a finite difference will be wrong by seven orders of magnitude for a forward difference. If you need machine precision, hand-code the derivative.

Measured. accuracy_floor_at_the_default_step in differentiate/tests.rs, release, 2026-08-13. Worst relative error over six points in [0.25, 3.3], each scheme at its own default step:

Function	forward	backward	central	central-4th
`sin`	1.281191e-8	1.401018e-8	6.528067e-11	1.785239e-13
`exp`	2.413003e-8	2.383524e-8	6.640831e-11	1.800352e-13
`x^3 - 2x`	4.470348e-8	4.470348e-8	4.583556e-11	9.992007e-14
`1/(1 + x^2)`	8.046627e-9	8.430242e-9	2.949979e-11	9.620083e-14
`tanh`	6.775837e-9	8.125324e-9	1.244682e-11	1.963985e-13
worst	4.470348e-8	4.470348e-8	6.640831e-11	1.963985e-13
predicted	1.490116e-8	1.490116e-8	3.666853e-11	3.000214e-13

Every scheme lands within a factor of three of its prediction, and the ordering is exactly as predicted: central-4th beats central by 338x, and central beats forward by 673x. The theory is a usable bound, not a story.

Observed convergence order. observed_convergence_order_matches_theory, same run. Absolute error of d/dx sin(x) at x = 1 (exact 5.40302305868139765e-1) against the relative step:

Relative step	forward	backward	central	central-4th
`1e-1`	4.293855e-2	4.113845e-2	9.000537e-4	1.125295e-7
`1e-2`	4.216325e-3	4.198315e-3	9.004993e-6	1.126843e-11
`1e-3`	4.208255e-4	4.206454e-4	9.005042e-8	1.049161e-13
`1e-4`	4.207445e-5	4.207265e-5	9.003700e-10	3.312906e-13
observed order	1.0079	0.9912	1.9998	3.9994

The fourth-order column turning back upward at 1e-4 (1.05e-13 to 3.31e-13) is the round-off wall arriving, exactly where eps^(1/5) = 7.4e-4 says it should.

The round-off wall. a_step_far_below_the_optimum_is_worse_not_better, same run. Central difference, d/dx sin(x) at x = 1, relative error:

Relative step	Relative error
`1.0000e-2`	1.666658e-5
`6.0555e-6` (`CBRT_EPSILON`)	5.373555e-12
`1.0000e-8`	5.303737e-9
`1.0000e-10`	1.909780e-7
`1.0000e-12`	5.609880e-5
`1.0000e-14`	7.666335e-3

A step six million times smaller than the optimum is ten million times worse: 1e-12 against 6.06e-6 moves the error from 5.373555e-12 to 5.609880e-5, a factor of 1.04e7. “Make h tiny for accuracy” is the intuition this table exists to destroy.

§Forward against central: the cost/accuracy trade, measured

For an n-dimensional Jacobian the evaluation counts are n + 1, 2n and 4n (see DiffScheme::evaluations_per_jacobian). The measured accuracy ratio from the table above is 673x for one extra evaluation per column going forward to central, and a further 338x for two more going to Central4th. Central is very nearly always the right default, which is why it is what DiffSettings::central exists for and what NumericalJacobian’s documentation recommends; DiffScheme::Forward is for the case where the function is genuinely expensive and 1e-8 is enough.

§Verification against analytic Jacobians

Methodology. Three systems whose Jacobians can be written down exactly are differenced and compared entry by entry: a quadratic system [x0^2 + x1, x0*x1^2], a trigonometric system [sin(x0)cos(x1), exp(x0)*x1] in which no derivative of any order vanishes, and a stiff linear pair [-1000*y0 + y1, y0 - y1] with a 1000:1 entry spread. Pass criterion: worst relative error below 1e-7 for the first-order schemes and 1e-10/1e-11 for the higher-order ones. All four schemes pass on all three systems (release, 2026-08-13).

The stiff pair is the informative one, because a linear system has exactly zero truncation error — so every digit lost is cancellation, and the measurement isolates it. Absolute error per entry at y = [0.4, -0.9]:

Scheme	`J[0][0]` (-1000)	`J[0][1]` (1)	`J[1][0]` (1)	`J[1][1]` (-1)
forward	0	0	0	0
backward	0	0	0	0
central	6.600658e-10	1.778424e-9	9.167112e-12	1.833422e-11
central-4th	5.456968e-12	6.380474e-11	2.498002e-14	0

J[0][1] is the worst entry by two orders of magnitude, and it is the small entry in the row that also holds -1000. Row 0 evaluates to about -400.9; perturbing x1 changes it by 1.2e-5, so forming that difference discards eleven significant digits before the division. This is a general property of finite-difference Jacobians, not of this implementation: an entry is only as accurate as its own magnitude relative to the largest term in its row. A badly-scaled system loses precision in exactly the entries a stiff solver most needs. If that matters, hand-code the Jacobian or rescale the equations.

The one-sided schemes returning exactly zero error here is a property of this particular linear case (f(x+h) - f(x) is exact when f is affine and the operands round identically) and must not be read as one-sided differences being more accurate in general — the accuracy table above shows them 673x worse on smooth non-linear functions.

§Measured cost against a hand-coded Jacobian

Methodology. Van der Pol mu = 5, y(0) = [2, 0], integrated over x in [0, 10] by crate::ode::Rosenbrock23 with abs_tol = rel_tol = 1e-8; best of 5, release, default features, the loaded 4-core machine described on DERIVATIVE_BATCH_MIN_POINTS. Van der Pol has an analytic Jacobian, which is the baseline. Produced by the #[ignore]d numerical_jacobian_overhead_benchmark, two runs.

Jacobian	Time (A)	Time (B)	vs analytic (A)	(B)	`y0(10)`
analytic (hand-coded)	849.99 us	855.47 us	1.00x	1.00x	-1.15870127
`DiffScheme::Forward`	1619.03 us	1603.88 us	1.90x	1.87x	-1.15870127
`DiffScheme::Backward`	1607.45 us	1614.33 us	1.89x	1.89x	-1.15870127
`DiffScheme::Central`	1673.66 us	1677.88 us	1.97x	1.96x	-1.15870127
`DiffScheme::Central4th`	2487.96 us	2455.24 us	2.93x	2.87x	-1.15870127

All four schemes reproduce the analytic result to all eight printed decimals, at roughly twice the cost. That is the honest headline: a numerical Jacobian is not free, it is not exact, and for a n = 2 stiff system it costs about a factor of two and changes nothing you can see.

§A NaN Jacobian is NOT reported by the solver — check the counter

When a Jacobian cannot be differenced, this module writes NaN into the entries and says so through DiffStatus. The natural expectation is that Rosenbrock23 then fails loudly. It does not, and this was measured rather than assumed (a_jacobian_that_cannot_be_differenced_is_counted_and_reaches_the_solver_as_nan):

crate::ode::Rosenbrock23::integrate returns Ok(());
the state vector comes back NaN.

The cause is in the ODE layer, not here: ode::normalize_error folds the per-equation errors with f64::max, and f64::max(0.0, NaN) is 0.0 — so a NaN error estimate looks like a perfectly converged step and every sub-step is accepted. Nothing in this module can change that.

Consequence for callers: NumericalJacobian::non_finite_jacobians is the only in-band signal that anything went wrong. Check it after any integration whose result you intend to trust. It is not decoration.

§Hybrid means dispatch, not two APIs

Every entry point takes a ComputeBackend parameter; there is no _parallel() sibling. With the parallel feature off, ComputeBackend::CpuMulti resolves down to ComputeBackend::Serial and the answer is unchanged — bit for bit, not merely close. There is no Gpu kernel here yet, so a Gpu request degrades to the best available CPU path.

Two independent parallel axes live here, with separately measured crossovers 2048x apart:

Entry point	Parallel over	Crossover
`derivative_batch`	independent points	`DERIVATIVE_BATCH_MIN_POINTS` = 65 536
`jacobian_batch`	independent lanes	`JACOBIAN_BATCH_MIN_PROBLEMS` = 256
`jacobian`	the columns of one Jacobian	`JACOBIAN_COLUMN_MIN_DIMENSION` = 32

They are not nested: jacobian_batch runs each lane’s columns serially, because the lane axis is already saturating the pool.

§Determinism — bitwise identical across backends and thread counts

This module returns bit-for-bit identical output on ComputeBackend::Serial and ComputeBackend::CpuMulti, at any thread count, on every run, provided the caller’s function is a deterministic pure function of its arguments.

The argument is the same one crate::math::minimise makes: lane i’s (or column j’s) answer is a pure function of its own samples, and no arithmetic crosses lanes or columns. A parallel sum would have to re-associate, and floating-point addition is not associative; a set of independent difference quotients has nothing to re-associate. Both backends call the same #[inline] per-lane and per-column kernels, and only the identity of the calling thread differs.

Verified by the bitwise_* tests in differentiate/tests.rs on 2 048 derivative lanes, 512 four-dimensional Jacobian lanes and one 96-dimensional Jacobian, all built with points spread over seven decades so the per-lane step differs. Measured 2026-08-13 (release, --features parallel, 4 logical cores): bit-identical on every observable field of every lane, for all four DiffScheme variants, at 1, 2, 4 and 8 workers. The single-point derivative and single-lane jacobian forms are separately asserted bit-identical to their one-element batches.

The #[ignore]d differentiate_thread_scaling_benchmark re-asserts the same identity while timing it, on 65 536 lanes with DiffScheme::Central4th (4 evaluations per lane), best of 7, two runs:

Worker threads	Time (A)	Speed-up (A)	(B)	Bitwise vs serial
serial reference	5902.59 us	1.00x	1.00x	—
1	6008.47 us	0.98x	0.98x	identical
2	3181.82 us	1.86x	1.85x	identical
4	1582.50 us	3.73x	3.70x	identical
8	1577.90 us	3.74x	3.83x	identical

The “identical” column is asserted by the benchmark, not merely printed. Scaling is close to linear to 4 workers and flat beyond, which is what four logical cores should do. The machine was not idle (see DERIVATIVE_BATCH_MIN_POINTS for the load); one machine, one batch, two runs, nothing measured on Android hardware or a many-core server.

The one way a caller can break this is to supply a function that is not pure — one that reads a random number generator, accumulates into shared interior-mutable state, or depends on the calling thread. The Sync bound permits it; this contract forbids it.

§Failure is reported, never swallowed

Every lane and every Jacobian carries a DiffStatus.
DerivativeSolution::derivative and JacobianSolution::matrix return Option, Some only on success. The diagnostic values are behind the deliberately-named DerivativeSolution::raw_value and JacobianSolution::raw_matrix.
DerivativeBatch::values and JacobianBatch::matrices are all-or-nothing: they return DiffBatchFailure naming the failure count and the first failing lane, rather than a Vec with a NaN in it.
A failed Jacobian column is NaN, never 0.0, and JacobianSolution::first_bad_column names the offending variable.

But read the limits of that guarantee — the “What is detected, and what cannot be” section on derivative lists three classes of bad input that no finite-difference kernel can detect, chief among them a singularity that a symmetric stencil steps over.

§What is deliberately NOT here: dual-number autodiff

Bead op-yvj.4.6 offers forward-mode dual numbers as an optional exact alternative, “only if it stays simple”. It is not implemented, on purpose. Making it useful means every function a caller wants differentiated must be generic over the scalar type, which would push a type parameter through OdeSystem, through the thermophysics kernels, and into every caller’s own code. That is precisely the rise in reader context load the crate-level “Human interface layer” rule forbids, and the bead itself ranks that rule above the convenience. A caller who wants exact derivatives should hand-code them; that is what crate::ode::OdeSystem::jacobian is for.

§Units

Everything here is dimensionless f64, and that is a deliberate decision rather than uom being stripped.

A derivative changes dimension. d(enthalpy)/d(temperature) is a heat capacity; d(pressure)/d(volume) is none of the three. A generic differentiator therefore has no single uom type it could return — the output type is a function of two input types, which Rust can only express through a trait with an associated output type, i.e. exactly the generic machinery the “Human interface layer” rule forbids adding for its own sake. The bead anticipates this and directs that, where the generic form cannot be typed cleanly, a small number of concrete typed wrappers is preferred over one generic nobody can read.

So: uom typing is applied at the boundary, by the caller — convert in, convert out — exactly as crate::math::minimise and crate::math::parallel do. The one place a dimension does appear in this module’s own API is DiffSettings::min_scale, which carries the units of the variable being perturbed; its documentation says so, because a caller differentiating with respect to a pressure in pascals near zero must not leave it at 1.0.

§Cargo features and portability

The rayon paths sit behind the crate’s parallel feature, which is off by default; with it off this module still compiles and every entry point still works. rayon is pure Rust with no system component, so everything here compiles and runs on aarch64-linux-android / Termux exactly as on desktop. Nothing in this module is target-gated.

§Example

use outram_foam_basic_lib::math::differentiate::{DiffSettings, NumericalJacobian};
use outram_foam_basic_lib::ode::{OdeSystem, Rosenbrock23};

// A stiff system with NO hand-coded Jacobian. Without the wrapper the
// default `OdeSystem::jacobian` would panic inside Rosenbrock23.
struct StiffPair;
impl OdeSystem for StiffPair {
    fn n_eqns(&self) -> usize { 2 }
    fn derivatives(&self, _x: f64, y: &[f64], dydx: &mut Vec<f64>) {
        dydx.clear();
        dydx.push(-1000.0 * y[0] + y[1]);
        dydx.push(y[0] - y[1]);
    }
}

let system = NumericalJacobian::new(StiffPair, DiffSettings::central());
let mut solver = Rosenbrock23::new(2, 1e-10, 1e-10);
let mut y = vec![1.0_f64, 1.0];
let mut dx = 1e-6;
solver.integrate(&system, 0.0, 1.0, &mut y, &mut dx).expect("integrates");

// The fast mode has decayed; the slow mode (eigenvalue about -0.999) remains.
assert!(y[1].abs() < 1.0 && y[1] > 0.0, "y1 = {}", y[1]);
// ALWAYS check this -- the solver does not report a NaN Jacobian itself.
assert_eq!(system.non_finite_jacobians(), 0);

Structs§

DerivativeBatch: The result of derivative_batch — one DerivativeSolution per point, in point order.
DerivativeSolution: One lane’s scalar derivative, with the diagnostics needed to judge it.
DiffBatchFailure: One or more lanes of a DerivativeBatch or JacobianBatch failed.
DiffSettings: Step-size policy for every entry point in this module.
JacobianBatch: The result of jacobian_batch — one JacobianSolution per lane.
JacobianSolution: One lane’s Jacobian, with the status needed to judge it.
NumericalJacobian: Wrap any OdeSystem so that crate::ode::Rosenbrock23 can integrate it without a hand-coded Jacobian.

Enums§

DiffScheme: Which finite-difference stencil to use.
DiffStatus: Why a derivative or Jacobian entry is, or is not, trustworthy.

Constants§

CBRT_EPSILON: f64::EPSILON.cbrt() = 6.0554544523933395e-6 — the relative step that balances truncation against round-off for a central difference.
DERIVATIVE_BATCH_MIN_POINTS: Point count below which a ComputeBackend::CpuMulti request runs derivative_batch on the calling thread instead.
FIFTH_ROOT_EPSILON: f64::EPSILON.powf(0.2) = 7.40095979741405e-4 — the relative step that balances truncation against round-off for the fourth-order Richardson scheme DiffScheme::Central4th.
JACOBIAN_BATCH_MIN_PROBLEMS: Lane count below which a ComputeBackend::CpuMulti request runs jacobian_batch on the calling thread instead.
JACOBIAN_COLUMN_MIN_DIMENSION: Dimension below which jacobian computes one Jacobian’s columns on the calling thread rather than spreading them across rayon.

Functions§

derivative: Differentiate one scalar function at one point.
derivative_backend_for: The ComputeBackend derivative_batch would actually use for n points if asked for requested — without running anything.
derivative_batch: Differentiate N independent scalar functions at N points, on the chosen backend.
jacobian: Assemble the Jacobian J[i][j] = d f_i / d x_j of one n-dimensional vector function at one point, by finite differences.
jacobian_batch: Assemble N independent Jacobians, one per lane, on the chosen backend.
jacobian_batch_backend_for: The ComputeBackend jacobian_batch would actually use for n independent Jacobian problems — without running anything.
jacobian_column_backend_for: The ComputeBackend jacobian would actually use to spread the columns of one dimension-dimensional Jacobian — without running anything.
ode_system_jacobian: Fill an OdeSystem’s Jacobian slots by finite differences.