Crate bulirsch 

Implementation of the Bulirsch-Stoer method for stepping ordinary differential equations.

The (Gragg-)Bulirsch-Stoer algorithm combines the (modified) midpoint method with Richardson extrapolation to accelerate convergence. It is an explicit method that does not require Jacobians.

This crate’s implementation contains simplistic adaptive step size routines with order estimation. Its API is designed to be useful in situations where an ODE is being integrated step by step with a prescribed time step, for example in integrated simulations of electromechanical control systems with a fixed control cycle period. Only time-independent ODEs are supported, but without loss of generality (since the state vector can be augmented with a time variable if needed).

The implementation follows:

• Press, William H. Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge University Press, 2007. Ch. 17.3.2.
• Deuflhard, Peter. “Order and stepsize control in extrapolation methods.” Numerische Mathematik 41 (1983): 399-422.

As an example, consider a simple oscillator system:

// Define ODE.
struct OscillatorSystem {
    omega: f64,
}

impl bulirsch::System for OscillatorSystem {
    type Float = f64;

    fn system(
        &self,
        y: bulirsch::ArrayView1<Self::Float>,
        mut dydt: bulirsch::ArrayViewMut1<Self::Float>,
    ) {
        dydt[[0]] = y[[1]];
        dydt[[1]] = -self.omega.powi(2) * y[[0]];
    }
}

let system = OscillatorSystem { omega: 1.2 };

// Set up the integrator.
let mut integrator = bulirsch::Integrator::default()
    .with_abs_tol(1e-8)
    .with_rel_tol(1e-8)
    .into_adaptive();

// Define initial conditions and provide solution storage.
let delta_t: f64 = 10.2;
let mut y = ndarray::array![1., 0.];
let mut y_next = ndarray::Array::zeros(y.raw_dim());

// Integrate for 10 steps.
let num_steps = 10;
for _ in 0..num_steps {
    integrator
        .step(&system, delta_t, y.view(), y_next.view_mut())
        .unwrap();
    y.assign(&y_next);
}

// Ensure result matches analytic solution.
approx::assert_relative_eq!(
    (system.omega * delta_t * num_steps as f64).cos(),
    y_next[[0]],
    epsilon = 5e-7,
    max_relative = 5e-7,
);

// Check integration performance.
assert_eq!(integrator.overall_stats().num_system_evals, 3770);
approx::assert_relative_eq!(integrator.step_size().unwrap(), 2.10, epsilon = 1e-2);

Note that 3.7k system evaluations have been used. By contrast, the ode_solvers::Dopri5 algorithm uses more:

struct OscillatorSystem {
    omega: f64,
}

impl ode_solvers::System<f64, ode_solvers::SVector<f64, 2>> for OscillatorSystem {
    fn system(
        &self,
        _x: f64,
        y: &ode_solvers::SVector<f64, 2>,
        dy: &mut ode_solvers::SVector<f64, 2>,
    ) {
        dy[0] = y[1];
        dy[1] = -self.omega.powi(2) * y[0];
    }
}

let omega = 1.2;
let delta_t: f64 = 10.2;
let mut num_system_eval = 0;
let mut y = ode_solvers::Vector2::new(1., 0.);
let num_steps = 10;
for _ in 0..num_steps {
    let system = OscillatorSystem { omega };
    let mut solver = ode_solvers::Dopri5::new(
        system,
        0.,
        delta_t,
        delta_t,
        y,
        1e-8,
        1e-8,
    );
    num_system_eval += solver.integrate().unwrap().num_eval;
    y = *solver.y_out().get(1).unwrap();
}
assert_eq!(num_system_eval, 7476);

// Ensure result matches analytic solution.
approx::assert_relative_eq!(
    (omega * delta_t * num_steps as f64).cos(),
    y[0],
    epsilon = 5e-7,
    max_relative = 5e-7,
);

As of writing this, the latest version of ode_solvers, 0.6.1, also gives a dramatically incorrect answer likely due to a regression. As a result we use version 0.5 as a dev dependency.

Structs

AdaptiveIntegrator

An ODE integrator using the Bulirsch-Stoer algorithm with an adaptive step size and adaptive order.

Integrator

An ODE integrator using the Bulirsch-Stoer algorithm with a fixed step size.

Stats

Statistics from taking an integration step.

StepSizeUnderflow

Error generated when integration produced a step size smaller than the minimum allowed step size.

Traits

Float

System

Trait for defining an ordinary differential equation system.

Type Aliases

ArrayView1

one-dimensional array view

ArrayViewMut1

one-dimensional read-write array view