Skip to main content

derivative

Function derivative 

Source
pub fn derivative<F>(x: f64, settings: DiffSettings, f: F) -> DerivativeSolution
where F: Fn(f64) -> f64,
Expand description

Differentiate one scalar function at one point.

The single-lane form of derivative_batch, for callers with one derivative to take. It runs on the calling thread — there is nothing to spread — and calls the same per-lane kernel, so it agrees with a one-element batch bit for bit.

§The realised-step correction

x + h is generally not representable, so the value the machine actually evaluates at differs from x + h in the last bits and the true step is not h. This kernel therefore evaluates at xp = x + h and divides by xp - x, which is exact, rather than by h. The device is from Numerical Recipes (Press et al., 3rd ed., section 5.7) and it removes an error source that would otherwise be comparable to the round-off term the step rule is trying to balance. DerivativeSolution::realised_step reports the corrected denominator.

§Arguments

  • x — the point, in the caller’s own units.
  • settings — scheme and step-size policy; see DiffSettings.
  • f — the function. Units of the return value are the caller’s.

§Returns

A DerivativeSolution whose derivative is Some only if every evaluation was finite and the realised step was non-zero. That is the whole of the guarantee — read the next section before relying on it.

§What is detected, and what cannot be

The status is computed from exactly one predicate — every sampled value and the resulting quotient are finite, and the realised step is non-zero — and DerivativeSolution::derivative returns Some on exactly that same predicate. They cannot disagree.

Detected: a non-finite sample (DiffStatus::NotFinite), a non-finite point (DiffStatus::InvalidPoint), and a step that rounds away so x + h == x (DiffStatus::DegenerateStep).

Not detected, and not detectable by any finite-difference kernel:

  • A singularity the stencil steps over. 1/x at x = 0 is sampled by the central stencil at +h and -h, both perfectly finite, so it returns 1/h^2 with DiffStatus::Ok. The kernel never evaluates at the pole and has no way to learn it is there. A one-sided scheme does see this particular case, because it evaluates at x itself — but it has the mirror-image blind spot on the other side.
  • Cancellation that leaves a finite number with no correct digits. The quotient is a perfectly ordinary f64; nothing about it says how many of its bits survived. This is what the step-size rule exists to bound, and why the module documents an accuracy floor rather than a guarantee.
  • A function that is not differentiable at x. |x| at 0 returns 0 from the central stencil, confidently.

If the function may be singular or kinked, bracket it away from the trouble or check the result against a second scheme; the status field will not do it for you.

§Example

use outram_foam_basic_lib::math::differentiate::{
    derivative, DiffSettings, DiffStatus,
};

// d/dx sin(x) at x = 1 is cos(1).
let s = derivative(1.0, DiffSettings::central(), |x: f64| x.sin());
let d = s.derivative().expect("finite everywhere");
assert!((d - 1.0_f64.cos()).abs() < 1e-10, "got {d}");

// A sample that comes back non-finite IS reported: the central stencil for
// sqrt at x = 0 evaluates at -h, which is NaN.
let bad = derivative(0.0, DiffSettings::central(), |x: f64| x.sqrt());
assert_eq!(bad.status(), DiffStatus::NotFinite);
assert!(bad.derivative().is_none());

// But a pole the stencil STEPS OVER is not, and cannot be -- see
// "What is detected, and what cannot be" above.
let undetected = derivative(0.0, DiffSettings::central(), |x: f64| 1.0 / x);
assert_eq!(undetected.status(), DiffStatus::Ok);
assert!(undetected.derivative().is_some());