Struct SolverSensi 

pub struct SolverSensi<UserData, F, FS, const N: usize, const N_SENSI: usize> { /* private fields */ }

The ODE solver with sensitivities.

Type Arguments

• F is the type of the right-hand side function

• FS is the type of the sensitivities right-hand side function

• UserData is the type of the supplementary arguments for the right-hand-side. If unused, should be ().

• N is the “problem size”, that is the dimension of the state space.

• N_SENSI is the number of sensitivities computed

Implementations

impl<UserData, F, FS, const N: usize, const N_SENSI: usize> Solver<UserData, F, FS, N, N_SENSI>

where F: Fn(Realtype, &[Realtype; N], &mut [Realtype; N], &UserData) -> RhsResult, FS: Fn(Realtype, &[Realtype; N], &[Realtype; N], [&[Realtype; N]; N_SENSI], [&mut [Realtype; N]; N_SENSI], &UserData) -> RhsResult,

pub fn new( method: LinearMultistepMethod, f: F, f_sens: FS, t0: Realtype, y0: &[Realtype; N], y_s0: &[[Realtype; N]; N_SENSI], rtol: Realtype, atol: AbsTolerance<N>, atol_sens: SensiAbsTolerance<N, N_SENSI>, user_data: UserData, ) -> Result<Self>

Creates a new solver.

Examples found in repository

examples/oscillator_sensi.rs (lines 34-45)

fn main() {
    let y0 = [0., 1.];
    //define the right-hand-side
    fn f(_t: Realtype, y: &[Realtype; 2], ydot: &mut [Realtype; 2], k: &Realtype) -> RhsResult {
        *ydot = [y[1], -y[0] * k];
        RhsResult::Ok
    }
    //define the sensitivity function for the right hand side
    fn fs(
        _t: Realtype,
        y: &[Realtype; 2],
        _ydot: &[Realtype; 2],
        ys: [&[Realtype; 2]; N_SENSI],
        ysdot: [&mut [Realtype; 2]; N_SENSI],
        k: &Realtype,
    ) -> RhsResult {
        // Mind that when indexing sensitivities, the first index
        // is the parameter index, and the second the state variable
        // index
        *ysdot[0] = [ys[0][1], -ys[0][0] * k];
        *ysdot[1] = [ys[1][1], -ys[1][0] * k];
        *ysdot[2] = [ys[2][1], -ys[2][0] * k - y[0]];
        RhsResult::Ok
    }

    const N_SENSI: usize = 3;

    // the sensitivities in order are d/dy0[0], d/dy0[1] and d/dk
    let ys0 = [[1., 0.], [0., 1.], [0., 0.]];

    //initialize the solver
    let mut solver = SolverSensi::new(
        LinearMultistepMethod::Adams,
        f,
        fs,
        0.,
        &y0,
        &ys0,
        1e-4,
        AbsTolerance::scalar(1e-4),
        SensiAbsTolerance::scalar([1e-4; N_SENSI]),
        1e-2,
    )
    .unwrap();
    //and solve
    let ts: Vec<_> = (1..100).collect();
    println!("0,{},{}", y0[0], y0[1]);
    for &t in &ts {
        let (_tret, &[x, xdot], [&[dy0_dy00, dy1_dy00], &[dy0_dy01, dy1_dy01], &[dy0_dk, dy1_dk]]) =
            solver.step(t as _, StepKind::Normal).unwrap();
        println!(
            "{},{},{},{},{},{},{},{},{}",
            t, x, xdot, dy0_dy00, dy1_dy00, dy0_dy01, dy1_dy01, dy0_dk, dy1_dk
        );
    }
}

pub fn step( &mut self, tout: Realtype, step_kind: StepKind, ) -> Result<(Realtype, &[Realtype; N], [&[Realtype; N]; N_SENSI])>

Takes a step according to step_kind (see StepKind).

Returns a tuple (t_out,&y(t_out),[&dy_dp(tout)]) where t_out is the time reached by the solver as dictated by step_kind, y(t_out) is an array of the state variables at that time, and the i-th dy_dp(tout) is an array of the sensitivities of all variables with respect to parameter i.

Examples found in repository

examples/oscillator_sensi.rs (line 52)

fn main() {
    let y0 = [0., 1.];
    //define the right-hand-side
    fn f(_t: Realtype, y: &[Realtype; 2], ydot: &mut [Realtype; 2], k: &Realtype) -> RhsResult {
        *ydot = [y[1], -y[0] * k];
        RhsResult::Ok
    }
    //define the sensitivity function for the right hand side
    fn fs(
        _t: Realtype,
        y: &[Realtype; 2],
        _ydot: &[Realtype; 2],
        ys: [&[Realtype; 2]; N_SENSI],
        ysdot: [&mut [Realtype; 2]; N_SENSI],
        k: &Realtype,
    ) -> RhsResult {
        // Mind that when indexing sensitivities, the first index
        // is the parameter index, and the second the state variable
        // index
        *ysdot[0] = [ys[0][1], -ys[0][0] * k];
        *ysdot[1] = [ys[1][1], -ys[1][0] * k];
        *ysdot[2] = [ys[2][1], -ys[2][0] * k - y[0]];
        RhsResult::Ok
    }

    const N_SENSI: usize = 3;

    // the sensitivities in order are d/dy0[0], d/dy0[1] and d/dk
    let ys0 = [[1., 0.], [0., 1.], [0., 0.]];

    //initialize the solver
    let mut solver = SolverSensi::new(
        LinearMultistepMethod::Adams,
        f,
        fs,
        0.,
        &y0,
        &ys0,
        1e-4,
        AbsTolerance::scalar(1e-4),
        SensiAbsTolerance::scalar([1e-4; N_SENSI]),
        1e-2,
    )
    .unwrap();
    //and solve
    let ts: Vec<_> = (1..100).collect();
    println!("0,{},{}", y0[0], y0[1]);
    for &t in &ts {
        let (_tret, &[x, xdot], [&[dy0_dy00, dy1_dy00], &[dy0_dy01, dy1_dy01], &[dy0_dk, dy1_dk]]) =
            solver.step(t as _, StepKind::Normal).unwrap();
        println!(
            "{},{},{},{},{},{},{},{},{}",
            t, x, xdot, dy0_dy00, dy1_dy00, dy0_dy01, dy1_dy01, dy0_dk, dy1_dk
        );
    }
}

Trait Implementations

impl<UserData, F, FS, const N: usize, const N_SENSI: usize> Drop for Solver<UserData, F, FS, N, N_SENSI>

fn drop(&mut self)

Executes the destructor for this type.