oscillator_sensi/
oscillator_sensi.rs

1use cvode_wrap::*;
2
3fn main() {
4    let y0 = [0., 1.];
5    //define the right-hand-side
6    fn f(_t: Realtype, y: &[Realtype; 2], ydot: &mut [Realtype; 2], k: &Realtype) -> RhsResult {
7        *ydot = [y[1], -y[0] * k];
8        RhsResult::Ok
9    }
10    //define the sensitivity function for the right hand side
11    fn fs(
12        _t: Realtype,
13        y: &[Realtype; 2],
14        _ydot: &[Realtype; 2],
15        ys: [&[Realtype; 2]; N_SENSI],
16        ysdot: [&mut [Realtype; 2]; N_SENSI],
17        k: &Realtype,
18    ) -> RhsResult {
19        // Mind that when indexing sensitivities, the first index
20        // is the parameter index, and the second the state variable
21        // index
22        *ysdot[0] = [ys[0][1], -ys[0][0] * k];
23        *ysdot[1] = [ys[1][1], -ys[1][0] * k];
24        *ysdot[2] = [ys[2][1], -ys[2][0] * k - y[0]];
25        RhsResult::Ok
26    }
27
28    const N_SENSI: usize = 3;
29
30    // the sensitivities in order are d/dy0[0], d/dy0[1] and d/dk
31    let ys0 = [[1., 0.], [0., 1.], [0., 0.]];
32
33    //initialize the solver
34    let mut solver = SolverSensi::new(
35        LinearMultistepMethod::Adams,
36        f,
37        fs,
38        0.,
39        &y0,
40        &ys0,
41        1e-4,
42        AbsTolerance::scalar(1e-4),
43        SensiAbsTolerance::scalar([1e-4; N_SENSI]),
44        1e-2,
45    )
46    .unwrap();
47    //and solve
48    let ts: Vec<_> = (1..100).collect();
49    println!("0,{},{}", y0[0], y0[1]);
50    for &t in &ts {
51        let (_tret, &[x, xdot], [&[dy0_dy00, dy1_dy00], &[dy0_dy01, dy1_dy01], &[dy0_dk, dy1_dk]]) =
52            solver.step(t as _, StepKind::Normal).unwrap();
53        println!(
54            "{},{},{},{},{},{},{},{},{}",
55            t, x, xdot, dy0_dy00, dy1_dy00, dy0_dy01, dy1_dy01, dy0_dk, dy1_dk
56        );
57    }
58}
oscillator_sensi/oscillator_sensi.rs

oscillator_sensi/
oscillator_sensi.rs
