Skip to main content

solve_forward_sensitivity

Function solve_forward_sensitivity 

Source
pub fn solve_forward_sensitivity<Sol, S, Sys>(
    system: &Sys,
    t0: S,
    tf: S,
    y0: &[S],
    options: &SolverOptions<S>,
) -> Result<SensitivityResult<S>, SolverError>
where Sol: Solver<S>, S: Scalar, Sys: ParametricOdeSystem<S>,
Expand description

Compute the forward sensitivity matrix S(t) = ∂y(t) / ∂p of an ODE solution using any Solver.

Builds the augmented state z = [y; vec(S)] (column-major sensitivity flattening; see module-level docs), drives Sol over the augmented system, and splits the trajectory back into y and sensitivity blocks at every output time the solver produced.

§Errors

Returns whatever error Sol::solve raises. If the underlying integration returns success = false (e.g. step-size collapse or max-steps reached), the returned SensitivityResult propagates that with empty trajectory vectors and the solver’s diagnostic message.

§Example

use numra_ode::sensitivity::{solve_forward_sensitivity, ParametricOdeSystem};
use numra_ode::{DoPri5, SolverOptions};

// dy/dt = -k * y, y(0) = 1, k = 0.5.
// Analytical: y(t) = exp(-k t), ∂y/∂k = -t exp(-k t).
struct Decay { k: f64 }
impl ParametricOdeSystem<f64> for Decay {
    fn n_states(&self) -> usize { 1 }
    fn n_params(&self) -> usize { 1 }
    fn params(&self) -> &[f64] { std::slice::from_ref(&self.k) }
    fn rhs_with_params(&self, _t: f64, y: &[f64], p: &[f64], dy: &mut [f64]) {
        dy[0] = -p[0] * y[0];
    }
    fn jacobian_y(&self, _t: f64, _y: &[f64], jy: &mut [f64]) { jy[0] = -self.k; }
    fn jacobian_p(&self, _t: f64, y: &[f64], jp: &mut [f64]) { jp[0] = -y[0]; }
    fn has_analytical_jacobian_y(&self) -> bool { true }
    fn has_analytical_jacobian_p(&self) -> bool { true }
}

let r = solve_forward_sensitivity::<DoPri5, _, _>(
    &Decay { k: 0.5 },
    0.0, 2.0,
    &[1.0],
    &SolverOptions::default().rtol(1e-8).atol(1e-10),
).unwrap();

let last = r.len() - 1;
let dy_dk = r.dyi_dpj(last, 0, 0);
let analytical = -r.t[last] * (-0.5 * r.t[last]).exp();
assert!((dy_dk - analytical).abs() < 1e-5);