Odeint library
| Mario Mulansky and Karsten Ahnert (2014), Scholarpedia, 9(12):32342. | doi:10.4249/scholarpedia.32342 | revision #145932 [link to/cite this article] |
The odeint (ordinary differential equation integration) library is a collection of advanced numerical algorithms to solve initial-value problems of ordinary differential equations. It is written in C++ using modern programming techniques to provide high generality at optimal performance.
Contents |
Introduction
The odeint library provides a number of algorithms to solve initial value problems of ODEs, as typically arising in Dynamical systems, Celestial mechanics, Computational Neuroscience and many other fields of science. The focus of the library are explicit methods, although some implicit routines are also available. Its main advantage over other ODE libraries is a strict separation of the numerical algorithms from the underlying arithmetical computations. This creates the possibility to use odeint on GPUs or HPC clusters, or to combine odeint with custom data types as well as other numerical libraries.
Odeint is written in C++ and uses modern programming techniques such as generic programming and template metaprogramming to ensure high flexibility at optimal performance. In 2012, odeint became part of the Boost C++ library collection. For this, it passed an extensive community review process, where the library's design, usability and documentation was positively evaluated by independent C++ experts. Odeint is actively maintained and expanded, and new releases are published as part of Boost several times a year.
Algorithms
The algorithms in odeint address initial value problems of first order ODEs defined in explicit form: \[ \dot{y} = f(t,y) , ~~ y(t_0) = y_0, \] where \(y\) is generally some \(N\) dimensional vector.
Odeint provides a number of numerical schemes to find an approximate solution for such ODEs. Those schemes can be divided into two main groups: explicit and implicit methods. The main focus of odeint lies on explicit methods, where several algorithm families are implemented, like explicit Runge-Kutta methods, extrapolation methods, linear multistep methods and others (see Table 1). For implicit methods, only (semi-) implicit Runge-Kutta schemes are provided.
Despite those mathematical families of algorithms, odeint introduces stepper categories corresponding to the functionality of the method. In particular, the following categories are used in odeint:
- Stepper
- Solvers of this category provide basic time stepping without step size control or dense output functionality. They simply iterate the solution of the ODE with constant step size.
- Controlled-Stepper
- These steppers are capable of performing step size control, which means adjusting the step-size to achieve a given accuracy of the numerical solution. This allows for larger step sizes in a controlled way and hence might speed up the numerical procedure. Controlled-Steppers are usually build upon an underlying Error-stepper.
- Dense-Output-Stepper
- Dense output steppers provide continuous interpolation between two consecutive time steps. The basic idea is to use a method with step size control to perform large steps and the use then interpolation to efficiently interpolate between two steps.
Table 1 gives an overview over the implemented routines, their mathematical family, category and properties.
| Method | Family | Category | Order | Step size control | Dense output |
|---|---|---|---|---|---|
| Euler | explicit Runge-Kutta | Stepper | 1 | no | no |
| Runge-Kutta4 | explicit Runge-Kutta | Stepper | 4 | no | no |
| Runge-Kutta-Cash-Karp | explicit Runge-Kutta | Controlled-Stepper | 5 | yes | no |
| Runge-Kutta-Dopri5 | explicit Runge-Kutta | Dense-Output-Stepper | 5 | yes | yes |
| Runge-Kutta-Fehlberg | explicit Runge-Kutta | Controlled-Stepper | 8 | yes | no |
| Adams-Bashforth | multistep | Stepper | configurable | no | no |
| Adams-Bashforth-Moulton | predictor-corrector multistep | Stepper | configurable | no | no |
| Bulirsch-Stoer | extrapolation | Dense-Output-Stepper | variable | yes | yes |
| Runge-Kutta-Nystroem | symplectic Runge-Kutta | Stepper | 4 | no | no |
| Implicit Euler | implicit Runge-Kutta | Stepper | 1 | no | no |
| Rosenbrock | semi-implicit Runge-Kutta | Dense-Output-Stepper | 4 | yes | yes |
| Velocit-Verlet | symplectic verlet | Stepper | 1 | no | no |
The Runge-Kutta4 stepper is the most renown algorithm, it provides a very robust and simple scheme to obtain an approximation using a fixed step-size. However, the more modern Adams-Bashforth-Moulton multi-step method provides a better efficiency with the same functionality and should be the preferred choice when a fixed step-size is reasonable. Often one can get even better performance by switching to a controlled stepper with dense-output like Runge-Kutta-Dopri5, where larger steps can be performed and intermediate results are interpolated in an efficient manner. If very high accuracy is required, the Bulirsch-Stoer method should be considered as it automatically adjusts the order as well as the step size.
Some situations require more specialized schemes, e.g. Runge-Kutta-Nystroem methods for Hamiltonian systems or the Velocity-Verlet scheme for molecular-dynamics and granular matter simulations.
Finally, for stiff systems one should employ implicit schemes such as the Rosenbrock method, which offers higher stability, but additionally requires the Jacobian \(Df\) of the rhs.
Features
Additionally to the basic solver algorithms listed above, odeint also provides some functionality helping to reduce the implementation effort for ODE simulations:
- integrate routines with observables and different observation strategies.
- iterator interface to the solvers.
- configurable data types.
- configurable computational backends.
- automatic memory management.
The fact that most algorithms in odeint have configurable data types and configurable computational backends ensures the library's high flexibility. For example, it is possible to work directly with ODEs defined for complex numbers, or to use arbitrary precision types if the usual double arithmetic is not sufficient. It is also possible to utilize the SIMD capabilities of modern CPUs, e.g. with help of the NT2 SIMD library within odeint. Furthermore, on can employ computations with physical units by using the Boost.Units library in combination with odeint. Another important point is its interoperability with linear algebra packages. In fact, odeint readily supports Eigen, Boost.Ublas, MTL4, and Blaze. Other frameworks can be adapted as well.
By virtue of exchanging the computational backends in odeint one can make use of parallelization capabilities. Odeint supports computational backends for the following parallel computing technologies:
- CUDA with the Thrust library
- OpenCL with the vexCL library
- OpenMP
- MPI
- HPX
Note, that for most situations the default configuration of odeint is sufficient and it is not necessary to select other backends.
Library Design
The odeint library is designed in a modularized way which allows to change certain parts of the implementation and to adjust the library to individual needs without re-implementing the fundamental algorithms.
The components of odeint are shown in Figure 1.
The core of the library are the Steppers that contain the actual algorithms for solving ODEs.
Additionally, odeint contains several integrate functions that iterate along a numerical solution of the given ODE.
The iterators provide a different interface for such iteration that can be used within standard algorithms, for example for_each, find or copy_if.
The computational backends of odeint are the Algebra and Operations.
They contain the functionality to perform the computations required by the algorithms in the Steppers.
Depending on the data structure used for representing the state of the ODE, different backends are requried.
A standard choice is to represent the state as a vector<double> which works with the default backend choice: range_algebra and default_operations.
However, the separation of stepper algorithm and computational backend allows to use the odeint library also in non-standard situations.
Odeint already includes several backends to work with the most common data structures and hardware setups.
Table 2 gives an exemplary (but not complete) overview of possible situation that can be readily addressed with the backends included in odeint.
| Data structure | Algebra | Operations | Remark |
|---|---|---|---|
std::vector, std::list, ...
|
range_algebra
|
default_operations
|
Standard situation: Containers |
boost::array
|
array_algebra
|
default_operations
|
System size known at compile time. |
double, complex
|
vector_space_algebra
|
default_operations
|
One-dimensional (complex) ODEs. |
ublas::vector, Eigen::Matrix, mtl::dense_vector,blaze::DynamicVector, ...
|
vector_space_algebra
|
default_operations
|
Linear algebra libraries. |
boost::fusion::vector
|
fusion_algebra
|
default_operations
|
Compile time sequence. |
thrust::vector
|
thrust_algebra
|
thrust_operations
|
GPU computation with CUDA. |
vexcl::vector
|
vector_space_algebra
|
default_operations
|
GPU computation with vexCL. |
odeint::openmp_state
|
openmp_algebra
|
default_operations
|
Parallel computations with OpenMP. |
odeint::mpi_state
|
mpi_algebra
|
default_operations
|
Parallel computations with MPI. |
Furthermore, with odeint, users can provide their own computational backend in case the problem requires some complex data structure not listed above.
An example could be ODEs on complex networks, which might require a complex_network_algebra implemented by the user.
Performance
Despite its modularized design and high flexibility, odeint still maintains competitive performance. Figure 2 shows a comparison of odeint with plain C and Fortran90 code. The performance is measured as the runtime required to perform 200,000 Runge-Kutta4 steps for the Lorenz system. On Intel hardware (Core i5 and Xeon E5), all implementations have virtually equivalent performance, while for the Opteron the Intel compiler produces the fastest code. Overall, this test shows that odeint provides competitive performance and that on modern compilers the abstraction and modularization induce no run-time costs.
On the contrary, the generic design of odeint makes it possible to achieve considerable performance gains, e.g. by utilizing the SIMD instructions of modern processors. For example, odeint readily works with NT2's Boost.SIMD library for which performance gains of up to a factor three on CPUs with AVX instruction sets have been observed for suitable problems.
Availability and Distribution
Odeint is distributed open-source under the Boost Licence. The odeint webpage contains all relevant information, source code links and documentation. New versions are published regularly as part of the Boost libraries several times a year. The latest version is available on Github, which also provides an issue system to submit bug reports.
References
- Karsten Ahnert and Mario Mulansky, odeint - Solving Ordinary Differential Equations in C++, AIP Conf. Proc. 1389 (ICNAAM 2011), pp. 1586-1589 (2011)
- Mario Mulansky and Karsten Ahnert, Generic Programming applied to Numerical Problems, AIP Conf. Proc. 1389 (ICNAAM 2011), pp. 1582-1585 (2011)
- Denis Demidov, Karsten Ahnert, Karl Rupp and Peter Gottschling, Programming CUDA and OpenCL: A Case Study Using Modern C++ Libraries. SIAM Journal on Scientific Computing, 35(5), C453-C472 (2013)
- Karsten Ahnert, Denis Demidov and Mario Mulansky, Solving Ordinary Differential Equations on GPUs, In Numerical Computations with GPUs, pp. 125-157. Springer International Publishing, 2014.
- Karsten Ahnert, odeint v2 - Solving ordinary differential equations in C++, Online Article (2011)
- Mario Mulansky, Boosting ODE simulations with Odeint and Boost.SIMD, Online Article (2014)
Internal References
- Kendall E. Atkinson (2007), Numerical analysis, Scholarpedia, 2(8):3163.
- Lawrence F. Shampine and Skip Thompson (2007), Initial value problems, Scholarpedia, 2(3):2861.
- John Butcher (2007), Runge-Kutta methods, Scholarpedia, 2(9):3147.
- Ernst Hairer and Gerhard Wanner (2010), Linear multistep method, Scholarpedia, 5(4):4591.
- Bill Gear (2007) Backward differentiation formulas, Scholarpedia, 2(8):3162.
- Lawrence F. Shampine and Skip Thompson (2007), Stiff systems, Scholarpedia, 2(3):2855.
- Philip Holmes and Eric T. Shea-Brown (2006) Stability, Scholarpedia, 1(10):1838.
- Alan C. Hindmarsh and Radu Serban (2007), Sundials equation solvers, Scholarpedia, 2(3):2860.
