ODEs Solving in Rust

Ode-solvers is a collection of numerical methods to solve ordinary differential equations (ODEs) in Rust. The following instructions should get you up and running in no time. See the examples and the API documentation for more details.

Installation

To start using the crate in your project, add the following dependency in your project's Cargo.toml file:

[dependencies]
ode_solvers = "0.5.x"

Then, in your main file, add

use ode_solvers::*;

Type alias definition

The numerical integration methods implemented in the crate support multi-dimensional systems. In order to define the dimension of the system, declare a type alias for the state vector. For instance

type State = Vector3<f64>;

The state representation of the system is based on the SVector<T, D> structure defined in the nalgebra crate. For convenience, ode-solvers re-exports six types to work with systems of dimension 1 to 6: Vector1<T>,..., Vector6<T>. For higher dimensions, the user should import the nalgebra crate and define a SVector<T, D> where the second type parameter of SVector is a dimension. For instance, for a 9-dimensional system, one would have:

type State = SVector<f64, 9>;

Alternativly, one can also use the DVector<T> structure from the nalgebra crate as the state representation. When using a DVector<T>, the number of rows in the DVector<T> defines the dimension of the system.

type State = DVector<f64>;

System definition

The system of first order ODEs must be defined in the system method of the System<T, V> trait. Typically, this trait is defined for a structure containing some parameters of the model. The signature of the System<T, V> trait is:

pub trait System<T, V> {
    fn system(&self, x: T, y: &V, dy: &mut V);
    fn solout(&self, _x: T, _y: &V, _dy: &V) -> bool {
        false
    }
}

where system must contain the ODEs: the second argument is the independent variable (usually time), the third one is a vector containing the dependent variable(s), and the fourth one contains the derivative(s) of y with respect to x. The method solout is called after each successful integration step and stops the integration whenever it is evaluated as true. The implementation of that method is optional. See the examples for implementation details.

Method selection

The following explicit Runge-Kutta methods are implemented in the current version of the crate:

Method Name Order Error estimate order Dense output order
Runge-Kutta 4 Rk4 4 N/A N/A
Dormand-Prince Dopri5 5 4 4
Dormand-Prince Dop853 8 (5, 3) 7

These methods are defined in the modules rk4, dopri5, and dop853. The Dormand-Prince methods feature:

  • Adaptive step size control
  • Automatic initial step size selection
  • Automatic stiffness detection
  • Sparse or dense output

The first step is to bring the desired module into scope:

use ode_solvers::dopri5::*;

Then, a structure is created using the new or the from_param method of the corresponding struct. Refer to the API documentation for a description of the input arguments.

let mut stepper = Dopri5::new(system, x0, x_end, dx, y0, rtol, atol);

The system is integrated using

let res = stepper.integrate();

which returns Result<Stats, IntegrationError>. Upon successful completion, res = Ok(Stats) where Stats is a structure containing some information on the integration process (such as number of function evaluations). If an error occurs, res = Err(IntegrationError). Finally, the results can be retrieved with

let x_out = stepper.x_out();
let y_out = stepper.y_out();

Adaptive Step Size

The adaptive step size used in Dopri5 and Dop853 is computed using a PI controller

where is a safety factor, is the current error , and is the previous error. The coefficients and determine the behavior of the controller. In order to make sure that doesn't change too abruptly, bounds on the factor of are enforced:

If , then ​ and if , .

The default values of these parameters are listed in the table below.

Dopri5 Dop853
0.17 1/8
0.04 0
0.2 0.333
10.0 6.0
0.9 0.9

These default values can be overridden by using the method from_param.

For Dopri5, a 5th order approximation ​and a 4th order approximation are constructed using the coefficients of the method. The error estimator is then

where is the dimension of the system and is a scale factor given by

For Dop853, an 8th order approximation , a 5th order approximation , and a 3rd order approximation are constructed using the coefficients of the method. Two error estimators are computed:

  and  

which are combined into a single error estimator

For a thorough discussion, see Hairer et al. [1].

References

[1] Hairer, E., Nørsett, S.P., Wanner, G., Solving Ordinary Differential Equations I, Nonstiff Problems, Second Revised Edition, Springer, 2008

[2] Hairer, E., Wanner, G., Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems, Second Revised Edition, Springer, 2002

[3] Butcher, J.C., Numerical Methods for Ordinary Differential Equations, Third Edition, John Wiley and Sons, 2016

[4] Gustafsson, K., Control-Theoretic Techniques for Stepsize Selection in Explicit Runge-Kutta Methods. ACM Transactions on Mathematical Software, 17, 4 (December 1991), 533-554.

[5] Gustafsson, K., Control-Theoretic Techniques for Stepsize Selection in Implicit Runge-Kutta Methods. ACM Transactions on Mathematical Software, 20, 4 (December 1994), 496-517.

[6] Söderlind, G., Automatic Control and Adaptive Time-Stepping, Numerical Algorithms (2002) 31: 281-310.