Crate russell_ode
source ·Expand description
Russell - Rust Scientific Library
russell_ode: Solvers for ordinary differential equations and differential algebraic equations
Important: This crate depends on external libraries (non-Rust). Thus, please check the Installation Instructions on the GitHub Repository.
The figure above is the solution of the Arenstorf Orbit obtained with russell_ode.
§Introduction
This library implements (natively) numerical solver solvers for systems of ordinary equations (ODEs) and differential-algebraic equation systems (DAEs) of Index-1. One advantage of a native implementation is the “safety aspects” enforced by Rust. Moreover, we implement thread-safe code. For example, the performance is improved when the real-based linear and complex-based linear systems are factorized concurrently, as in our Radau5.
The principal structs are:
- System defines the ODE or DAE system
- OdeSolver implements the “time-stepping” loop and calls the actual numerical solver
- Params holds numeric parameters needed by all methods
- (optional) Output holds the results from accepted steps (all methods) or the dense output (DoPri5, DoPri8, and Radau5 only)
- (optional) Stats holds statistics and benchmarking data
The System struct holds the number of equations (system dimension), the right-hand side function f(x, y), an optional function to compute the Jacobian matrix, and, also optionally, the mass matrix.
The OdeSolver approximates the solution of the ODE/DAE using either fixed or variable steps. Some methods can only be run with fixed steps (this is automatically detected). In addition to the system struct, the solver takes Params as input.
A set of default (~optimal) parameters are allocated by Params::new(). If needed, the user may tweak the parameters by accessing each parameter subgroup:
- ParamsNewton parameters for Newton’s iterations’ for the methods that use iterations such asBwEuler and Radau5
- ParamsStep parameters for the variable-step control
- ParamsStiffness parameters to control and enable the stiffness detection (DoPri5 and DoPri8 only)
- ParamsBwEuler parameters for the BwEuler solver
- ParamsRadau5 parameters for the Radau5 solver
- ParamsERK parameters for all explicit Runge-Kutta methods (e.g., DoPri5, DoPri8)
The ODE and DAE systems are represented as follows:
d{y}
[M] ———— = {f}(x, {y})
dx
where x is the independent scalar variable (e.g., time), {y} is the solution vector,
{f} is the right-hand side vector, and [M] is the so-called “mass matrix”.
Note: The mass matrix is optional and need not be specified (unless the DAE under study requires it).
The Jacobian is defined by:
∂{f}
[J](x, {y}) = ————
∂{y}
where [J] is the Jacobian matrix.
Note: The Jacobian function is not required for explicit Runge-Kutta methods (see Method and Information). Thus, one may simply pass the no_jacobian function and set HasJacobian::No in the system.
The flag ParamsNewton::use_numerical_jacobian may be set to true to compute the Jacobian matrix numerically. This option works with or without specifying the analytical Jacobian function.
§Recommended methods
- Method::DoPri5 for ODE systems and non-stiff problems using moderate tolerances
- Method::DoPri8 for ODE systems and non-stiff problems using strict tolerances
- Method::Radau5 for ODE and DAE systems, possibly stiff, with moderate to strict tolerances
Note: A Stiff problem arises due to a combination of conditions, such as the ODE system equations, the initial values, the stepsize, and the numerical method.
§Limitations
- Currently, the only method that can solve DAE systems is Method::Radau5
- Currently, dense output is only available for Method::DoPri5, Method::DoPri8, and Method::Radau5
§References
- E. Hairer, S. P. Nørsett, G. Wanner (2008) Solving Ordinary Differential Equations I. Non-stiff Problems. Second Revised Edition. Corrected 3rd printing 2008. Springer Series in Computational Mathematics, 528p
- E. Hairer, G. Wanner (2002) Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems. Second Revised Edition. Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
§Examples
See also the examples on GitHub page
§Simple system with mass matrix
Solve with Radau5:
y0' + y1' = -y0 + y1
y0' - y1' = y0 + y1
y2' = 1/(1 + x)
y0(0) = 1, y1(0) = 0, y2(0) = 0
Thus:
M y' = f(x, y)
with:
┌ ┐ ┌ ┐
│ 1 1 0 │ │ -y0 + y1 │
M = │ 1 -1 0 │ f = │ y0 + y1 │
│ 0 0 1 │ │ 1/(1 + x) │
└ ┘ └ ┘
The Jacobian matrix is:
┌ ┐
df │ -1 1 0 │
J = —— = │ 1 1 0 │
dy │ 0 0 0 │
└ ┘
Code:
use russell_lab::{vec_max_abs_diff, StrError, Vector};
use russell_ode::prelude::*;
use russell_sparse::CooMatrix;
fn main() -> Result<(), StrError> {
// DAE system
let ndim = 3;
let jac_nnz = 4; // number of non-zero values in the Jacobian
let mut system = System::new(
ndim,
|f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| {
f[0] = -y[0] + y[1];
f[1] = y[0] + y[1];
f[2] = 1.0 / (1.0 + x);
Ok(())
},
move |jj: &mut CooMatrix, alpha: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| {
jj.reset();
jj.put(0, 0, alpha * (-1.0)).unwrap();
jj.put(0, 1, alpha * (1.0)).unwrap();
jj.put(1, 0, alpha * (1.0)).unwrap();
jj.put(1, 1, alpha * (1.0)).unwrap();
Ok(())
},
HasJacobian::Yes,
Some(jac_nnz),
None,
);
// mass matrix
let mass_nnz = 5; // number of non-zero values in the mass matrix
system.init_mass_matrix(mass_nnz).unwrap();
system.mass_put(0, 0, 1.0).unwrap();
system.mass_put(0, 1, 1.0).unwrap();
system.mass_put(1, 0, 1.0).unwrap();
system.mass_put(1, 1, -1.0).unwrap();
system.mass_put(2, 2, 1.0).unwrap();
// solver
let params = Params::new(Method::Radau5);
let mut solver = OdeSolver::new(params, &system)?;
// initial values
let x = 0.0;
let mut y = Vector::from(&[1.0, 0.0, 0.0]);
// solve from x = 0 to x = 2
let x1 = 2.0;
let mut args = 0;
solver.solve(&mut y, x, x1, None, None, &mut args)?;
println!("y =\n{}", y);
Ok(())
}Modules§
- Makes available common structures and functions to perform computations
Structs§
- Holds information about the numerical method to solve (approximate) ODEs
- Implements a numerical solver for systems of ODEs
- Holds a counter of how many output files have been written
- Holds the data generated at an accepted step or during the dense output
- Holds the results at accepted steps or interpolated in a “dense” sequence of steps (dense output)
- Holds all parameters for the ODE Solver
- Holds the parameters for the BwEuler method
- Holds the parameters for explicit Runge-Kutta methods
- Holds parameters for the Newton-Raphson method
- Holds the parameters for the Radau5 method
- Holds parameters to control the variable stepsize algorithm
- Holds parameters for the stiffness detection algorithm
- Implements the Finite Difference (FDM) Laplacian operator in 2D
- Holds a collection of sample ODE problems
- Holds statistics and benchmarking data
- Defines a system of first order ordinary differential equations (ODE) or a differential-algebraic equations (DAE) of Index-1
Enums§
- Indicates whether the analytical Jacobian is available or not
- Specifies the numerical method to solve (approximate) ODEs
- Specifies the (boundary) side of a rectangle
Constants§
- Default number of steps to use when the automatic stepping is not available
Functions§
- Implements a placeholder function for when the analytical Jacobian is unavailable
Type Aliases§
- Defines a function of space that returns f64 (e.g., to calculate boundary condition values)
- Indicates that the system functions do not require extra arguments
- Defines a callback function to be executed during the output of results (accepted step or dense output)
- Defines the error output as a static string
- Defines a function to compute y(x) (e.g, the analytical solution)