## cvode-wrap

A wrapper around cvode and cvodeS from sundials, allowing to solve ordinary differential equations (ODEs) with or without their sensitivities

### 4 releases

 0.1.3 Jun 14, 2021 Jun 10, 2021 Jun 10, 2021 Jun 10, 2021

#317 in Science

BSD-3-Clause

33KB
656 lines

A wrapper around the cvode(S) ODE solver from sundials.

# Building sundials

To build sundials, activate the `sundials-sys/build_libraries` feature.

# Examples

## Oscillator

An oscillatory system defined by `x'' = -k * x`.

### Without sensitivities

``````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
}
//initialize the solver
let mut solver = SolverNoSensi::new(
f,
0.,
&y0,
1e-4,
AbsTolerance::scalar(1e-4),
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]) = solver.step(t as _, StepKind::Normal).unwrap();
println!("{},{},{}", t, x, xdot);
}
``````

### With sensitivities

The sensitivities are computed with respect to `x(0)`, `x'(0)` and `k`.

``````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(
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
);
}
``````

~9.5MB
~225K SLoC