Struct NumericalJacobian 

pub struct NumericalJacobian<S> { /* private fields */ }

Wrap any OdeSystem so that crate::ode::Rosenbrock23 can integrate it without a hand-coded Jacobian.

The problem this solves

OdeSystem::jacobian has a default body that is unimplemented!(), so a system that does not override it panics the moment a stiff solver asks for a Jacobian — inside Rosenbrock23::inner_step, and, if the integration is running in an ensemble, out through the rayon scope. Every system that only knows its own derivatives is locked out of the crate’s only stiff solver.

Wrapping it in NumericalJacobian supplies the missing method by finite differences and changes nothing else: n_eqns and derivatives are forwarded verbatim.

Owning, not borrowing

The wrapper owns the system by value, so it needs no lifetime parameter and no Box — both forbidden by the workspace design rules. Construct it with Self::new, get the system back with Self::into_inner.

An analytic Jacobian is still better

Finite differences cost n + 1 to 4n extra derivatives calls per Rosenbrock stage and are accurate to roughly sqrt(eps) to eps^(4/5) rather than to machine precision — see the module-level “Achievable accuracy” table. If the analytic Jacobian is available, write it. This wrapper is for the systems where it is not, and as a verification oracle for the ones where it is: differencing a system that also implements jacobian analytically and comparing is the cheapest real check that the hand-derived version has no sign or transposition error.

Units

Inherited from the wrapped system; nothing here is dimensioned.

Example

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

// A stiff scalar system with NO hand-coded Jacobian: dy/dx = -1000 y.
struct StiffDecay;
impl OdeSystem for StiffDecay {
    fn n_eqns(&self) -> usize { 1 }
    fn derivatives(&self, _x: f64, y: &[f64], dydx: &mut Vec<f64>) {
        dydx.clear();
        dydx.push(-1000.0 * y[0]);
    }
    // no `jacobian` override -- the default would panic
}

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

// exp(-1000 * 0.01) = exp(-10) = 4.5399929762484854e-5
let exact = (-10.0_f64).exp();
assert!((y[0] - exact).abs() < 1e-8, "got {}, want {exact}", y[0]);
assert_eq!(system.non_finite_jacobians(), 0);

Implementations

impl<S> NumericalJacobian<S>

pub fn new(system: S, settings: DiffSettings) -> Self

Wrap system, differencing its derivatives with settings.

Arguments

• system — owned by value.
• settings — scheme and step-size policy. DiffSettings::central is the usual choice for a Rosenbrock Jacobian: 2n evaluations for O(h^2) truncation, where DiffScheme::Central4th doubles the cost again for accuracy the step controller cannot exploit.

pub fn inner(&self) -> &S

Borrow the wrapped system.

pub fn into_inner(self) -> S

Unwrap, returning the system.

pub fn settings(&self) -> DiffSettings

The step-size policy in force.

pub fn non_finite_jacobians(&self) -> usize

How many OdeSystem::jacobian calls have failed since construction.

The trait method returns (), so it has nowhere to report a failure; this counter is that report. A non-zero value means at least one Jacobian was handed to the solver with NaN entries, which the solver will have turned into crate::ode::OdeError::StepSizeUnderflow rather than a wrong answer — but it is worth knowing why an integration failed, and “the Jacobian could not be differenced” is a different bug from “the system is too stiff”.

Units

A count, dimensionless.

Trait Implementations

impl<S: Debug> Debug for NumericalJacobian<S>

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl<S> OdeSystem for NumericalJacobian<S>

where S: OdeSystem,

fn n_eqns(&self) -> usize

Number of coupled equations (the length of the state vector y).

fn derivatives(&self, x: f64, y: &[f64], dydx: &mut Vec<f64>)

Fill dydx with the derivatives at (x, y).

fn jacobian( &self, x: f64, y: &[f64], dfdx: &mut Vec<f64>, dfdy: &mut SquareMatrix, )

Fill dfdx and dfdy with the Jacobian at (x, y).