pub fn solve_forward_sensitivity<Sol, S, Sys>(
system: &Sys,
t0: S,
tf: S,
y0: &[S],
options: &SolverOptions<S>,
) -> Result<SensitivityResult<S>, SolverError>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);