Introduction

This article introduces the second version of odeint - a C++ framework for solving ordinary differential equation (ODEs). It is based on template metaprogramming, is independent of a specific container type and can be used with modern graphic cards. Under to hood of odeint a lot of stuff has changed in the current version. It is not directly compatible with the first version but only slight changes have to be done to make version one and version two compatible. A brief introduction of the first version of odeint has been given here at codeproject.

The current version of odeint-v2 can be obtained from github. The full documentation includes a lot of tutorials and example and also can be found at github. odeint has been supported by the Google summer of code program. It is planned to release the odeint within the boost libraries, although it is not clear if odeint will pass the review process of boost.

Background

Solving ordinary differential equations is a very import task in mathematical modeling of physical, chemical, biological and even social systems. There is a long tradition of analyzing the methods of solving ODEs. These methods are summarized within the mathematical discipline Numerical Analysis. An ordinary differential equation has always the form:

d x(t) / d t = f( x(t) , t )

x is said to be the state of the ODE and can be vector wise. t is the dependent variable. We will call t always time although it can have a different (physical) meaning. The function f defines the ODE. Associated with every ODE is an initial value problem (IVP) that is given by the ODE plus an initial value x(t0)=x0. The solution of the initial value problem is the temporal evolution of x(t), with the additional condition that x(t0)=x0. It can be shown that for most of the physically relevant ODEs, every IVP has a unique solution. For more details see the background section in the article of the first version or have a look at some classical textbooks about differential equations or numerical analysis.

What's new?

Introduction of algebras and operations. They let you change the way how the basic mathematical operations are performed. With these concepts it is possible to solve ordinary differential equations with CUDA or by using sophisticated numerical libraries like the Intel Math Kernel library or the SIMD library.

Dense output. Methods with dense output are now supported by odeint. Dense output means that the solution is interpolated between intermediate points usually with the same order as the numerical solution. An own concept for dense output steppers has been introduced and the integrate functions support methods with dense output.

Meta functions and classes for the classical Runge-Kutta steppers. A new concept of generalizing the Runge-Kutta solvers has been introduced. This concept is based heavily on template metaprogramming. A specific solver is only defined by its coefficients. The performance of this methods is equal to the hand-written versions.

New methods. New methods has been introduced, like the Dormand-Prince 5(4) method, solvers for stiff systems like Rosenbrock-4 and multi step methods like Adams-Bashforth-Moulton. One attempt also tries to implement the Taylor series method.

Support for fancy libraries. odeint plays well with Boost.Range and Boost.Units. The first one lets you solve only a part of a full state while the second library performs compile-time dimension checks of a mathematical equation.

Factory functions. For easy generation of the controlled steppers and dense-output steppers.

Examples

Before going into the details of odeint we will show some examples.

Integration of the Lorenz system

Lets start with a examples from the previous articles -- the famous Lorenz system with its butterfly-like wing structure. Using odeint is very simple, just include the appropriate header files. odeint is completely header only. You will not have to link against other libs. This has the advantage that the compiler can do powerful optimization techniques resulting in a very fast run-time code. To use odeint just include

#include<boost/numeric/odeint.hpp>

The whole library lives in the namespace boost::numeric::odeint :

usingnamespace boost::numeric::odeint;

The definition of the Lorenz system is again either a function or a functor. But before declaring the ODE we need to define the state type, hence the type which odeint uses inside its steppers to store intermediate values and which is used in the system functions.

Now, we can in principle start to solve the Lorenz system. First, we have to choose a stepper. Lets start with classical Runge Kutta stepper of fourth order. We integrate the system for 1000 steps with a step size of 0.01.

Besides the integrate_const function you can also call the stepper directly. This might have some advantages if you want to perform non-trivial code between consecutive steps, but we recommend you to use the integrate functions and observers whenever possible. In the case of controlled and dense output steppers the integrate functions implement some complicate logic needed for a proper integration. Nevertheless the above code is exactly equivalent to

In most situation you want to do something with the state during the integration, for example write the state into a file or perform some fancy statistical methods. You can simply create an observer and call the integrate function with this observer.

In the previous examples the stepper was always a fourth order Runge Kutta stepper. This stepper is really nice. It is well known and widely used. But it does not have step size control not to mention dense output. If you want to use a controlled stepper which does an adaptive integration of you ODE you need a controlled stepper. For example you can use the Runge Kutta Cash Karp method

The runge_kutta_cash_karp54< state_type > is a simple error stepper which estimates the error made during one step. This stepper is used within the controlled_error_stepper which implements step size control on the basis of the error estimated from the Runge-Kutta Cash Karp stepper. integrate_adaptive performs the adaptive integration. This results in changes of the step size during time evolution.

Another error stepper is the so called Dormand-Prince 5(4) stepper that also possess dense-output functionality. To use it you need to nest the controller in a general dense output stepper:

Now, we have again used integrate_const, which uses all features of the dense output stepper. In principle one could also use a controlled stepper like that one from the previous example in integrate_const. But in this case the adaptive integration is performed exactly to the time points 0, dt, 2dt, 3dt, ... For small dt one would loose a lot of performance.

Of course, one can also tune the precision for the step size control, hence a measure which determines if the step made is to large or to small. For the example of the Runge-Kutta Cash-Karp stepper the precision is tuned via

abs_error determines the tolerance of the absolute error while rel_error is the tolerance of the relative error.

Factory functions

To ease the creation of the controlled steppers or the dense output steppers factory functions have been introduced. For example, a dense output stepper of the Dormand-Prince 5(4) stepper can easily be created by using:

The two parameters of the make_dense_output are the absolute and the relative error tolerances. For the controlled steppers the factory function is called make_controlled and can be used in exactly the same way as make_dense_output.

Solving ordinary differential equations on the GPU

The second example shows how odeint can be used to solve ordinary differential equations on the GPU. Of course not every problem is well suited to be ported on the GPU. For example the Lorenz in the last section is not a good candidate, it is simply to small. Using the GPU makes only sense for very large systems. One call on the GPU is slow, but this single call can perform a huge number of operations in parallel. Examples for problems which might suitable for the GPU are large lattices, like discretizations of partial differential equations (PDEs) or ensembles of ODEs.

In the example we will show how one can integrate a whole ensemble of uncoupled Lorenz systems where one parameter is varied. This example can be used for parameter studies. In the tutorial of the official documentation the same example is explain in more detail, here we only show how the GPU can be used.

At the moment odeint supports CUDA via Thrust. In principle it is also possible to use OpenCL, but this method is not yet implemented. Thrust is a STL like interface for the native CUDA API. If follows a functional programming paradigm. The standard way of using it is to call a general algorithm like thrust::for_each with a specific functor performing the basic operation. An example is

which generically adds 2 to every element of the container x. Another very nice feature of thrust is that the code can be used without modification with OpenMP, hence you can distribute your computations on the cores of your CPU. All you have to do to use this feature is to compile your code with -Xcompiler -fopenmp -DTHRUST_DEVICE_BACKEND=THRUST_DEVICE_BACKEND_OMP.

.

To use Thrust with odeint you have to:

use thrust::device_vector< double > as state type,

use a specialized algebra and specialized operations, namely thrust_algebra, thrust_operations (wait a second to see how this works),

Making your system function thrust-ready.

The first to points can easily be solved by defining the right state type and the right stepper:

The third point is the most difficult one, although thrust make this task fairly easy. We want to integrate an ensemble of N Lorenz systems, where every subsystem has a different parameter R. The dimension of a single Lorenz system is three, hence the dimension of the state type is 3*N. We will use a memory layout, where the first N entries in the state vector are the x components, the second N entries are the y components and the third N entries are the z components. Furthermore, we want to vary the parameter beta in the Lorenz system. Therefore, we also need a vector with N entries for the parameter. The layout of our Lorenz class is

Now, we are ready to compute dxdt. To do so we use a very nice feature of thrust - the zip iterator. Zip iterators allows us to iterate a whole bunch of vectors at once and perform one operation on this bunch.

So, we perform an iteration over x,y,z,R,dxdt,dydt,dzdt. The values are put inside a thrust::tuple. During the iteration lorenz_functor() is called with the tuple x,y,z,R,dxdt,dydt,dzdt. Everything we need to know about one single system is given in this tuple and we can calculate the the r.h.s. of the Lorenz system:

Of course for a real world application this example is not very interesting since nothing happens with this ensemble of Lorenz systems. In the tutorial section of odeint the same example is used to calculate the largest Lyapunov exponent for every single system. With the Lyapunov exponent it is possible to identify different regions in parameter space. For the Lorenz system regions with chaotic, regular (oscillating) and vanishing dynamics can be identified. The method of using thrust and CUDA for parameter studies is very efficient since many systems are solved in parallel. One can easily gain the performance by a factor of 20 or even more in comparison with the CPU.

Design and structure

In the previous section we have seen two examples of how odeint can be used to solve ordinary differential equations. This section is devoted to the technical details of odeint. Furthermore it gives an overview over existing methods and steppers in odeint.

odeint consists in principle of four parts:

Integrate functions

Steppers

Algebras and operations

Utilities

The integrate functions are responsible to perform the iteration of the ODE. They do the step size control and they might make use of dense output. The integrate functions come along with an observer concept which lets you observe your solution during the iteration. The steppers are the main building blocks of odeint. Here, the methods are implemented. Several steppers for different purposes exists, like

the classical Runge-Kutta steppers,

symplectic steppers for Hamiltonian systems,

implicit methods, and

stiff solvers.

The algebras and operations build an abstraction layer which lets you change the way how odeint performs basic numerical manipulations like addition, multiplication, etc. In the utility part you will find functionality like resizing and copying of state types. In the following we will introduce the several parts in more detail.

Steppers concepts

A stepper in odeint performs one single step. Several different stepper types exist in odeint. The most simple one is described by the Stepper concept. It has only one method do_step which performs the step. An example is the classical Runge-Kutta stepper of fourth order.

runge_kutta4< state_type > rk4;
rk4.do_step( system , x , t , dt );

Most of the stepper have a set of template parameters which tune their behavior. In the Runge-Kutta 4 example above you have explicitly to state the state type. This type is also used to store intermediate values. For the same stepper a lot of other parameter exist which are called with some default values. For example you can explicitly state the value type which gives the numerical values which is used during computations, the derivative type, the time type, the algebra and the operations. Furthermore, a policy parameter for resizing exists. In most cases the default parameters are a good choice. For exotic applications like using Boost.Units or for exotic algebras you need to specify them explicitly. The documentation of odeint shows some examples.

The first parameter of the do_step method specifies the ODE. It is expected that the system is either a function or functor and must be a callable object with the signature

system( x , dxdt , t )

Another concept is the ErrorStepper. Its basic usage is very similar to the Stepper concept:

The main difference is that it estimates the error made during one step and stores this error in xerr. Another concept is the ControlledStepper concept, which tries to perform one step. An example for the ControlledStepper concept is the controlled_runge_kutta which has one template argument, specifying the underlying Runge-Kutta stepper. For example, the Runge-Kutta-Cask-Karp stepper can be used via:

A stepper which models the ControlledStepper concept has usually some error bounds which must be fulfilled. If try_step is called the stepper will perform one step with a given step size dt and checks the error made during this step. If the error is within the desired tolerances it accepts the step and possibly increases the step size for the next evaluation. In case the error bounds are reached the step is rejected and the step size is decreased.

The fourth stepper concept is the DenseOutputStepper. A dense output stepper has two basic method do_step and calc_state. The first one performs one step. The state of the ODE is managed internally and usually this step is made adaptive, hence it is performed by a controlled stepper. The second function calculates intermediate values of the solution. Usually, the precision of the intermediate values are of the same order as the solution. An example is the Dormand-Prince stepper, which can be used like

Integrate functions

The integrate function are responsible for integrating the ODE with a specified stepper. Hence, they do the work of calling the do_step or try_step methods for you. This is especially useful if your are using a controlled stepper or a dense output stepper, since in this case you need some logic of calling the stepping method. Of course, you can always call the stepping methods by hand. If you are using a classical stepper without step size control and dense output functionality this might even be simpler than calling an integrate method.

odeint defines four different integrate methods.

integrate_const

integrate_adaptive

integrate times

integrate_n_steps

integrate_const performs the integration with a constant step size. At each step an observer can be called. An example is

integrate_const is implemented to use all features of the stepper. If you call it with a dense output stepper it will perform the stepping and calculate the state with the help of the dense output functionality. The same holds for the controlled steppers. Here the step size is adapted at each step.

The other integrate functions are similar to the integrate_const function. For example, integrate_adaptive performs adaptive integration of the ODE. integrate_times expects a range of times where you want to obtain the solution and integrate_n_steps integrates exactly n steps. This methods only works for the Stepper concept.

Algebras and operations

The algebra and operations introduce a degree of freedom which allows you to change and to adjust the way how odeint performs the basic numerical computations. In general an algebra defines how you can iterate over a set of state types while the operations perform the numerical operations like addition, subtraction, multiplication, etc. with the entries of the state type. Of course, the algebras and operations which can be used depend strongly on the state type.

Not all steppers provide the possibility to change algebras and operations. For example the Rosenbrock stepper and the implicit Euler stepper don't. They need to do matrix computations and thus depend on boost::numeric::ublas::vector and boost::numeric::ublas::matrix which provide basic vector and matrix algorithms. But the Runge-Kutta stepper as well as the symplectic and the multi-step methods use the algebras and operations.

Each algebra consists of a set of member functionsfor_each1( c1 , op ), for_each2( c1 , c2 , op ), for_each3( c1 , c2 , c3 , op ), ... and accumulate( c1 , op ). The variables c1, c2, c3, ... represent the state of the ODE or its first derivative. They are usually of the same type as the state_type of the stepper, but they must not. op is a function or functor which is performed on the elements of the state. In most cases this functor is taken from the operations. An operation is said to be a class with appropriate member structs scale_sum1, scale_sum2, scale_sum3, ... Each of this structs possesses an operator() that performs the operation. For more details about algebras and operations consult the documentation of odeint or read the source code.

The algebras and operations enter the stepper as additional template parameters. For example, all of the Runge-Kutta steppers have the following parameters:

Here, state_type represents the state of the ODE and value_type is the basic numerical type like double or complex. deriv_type represents the derivative of the state type hence the left hand side of the main equations while time_type is the type of the dependent variable. In most cases deriv_type is the same type as state_type and time_type as value_type. The next two parameters are then the algebra and the operations. All types except the state_type have default values chosen in a way suitable for most cases.

An example where all of these types are different is Boost.Units which can be used to implement compile-time validation of the dimension of an equations, see the tutorial in the odeint docs. Another example for algebras and operations different then the default one is the above example of thrust, which requires a special thrust_algebra as well as special thrust_operations. There are other algebras for example for the Intel math kernel library or you can easily implement your own algebra working with your own container.

The Runge-Kutta meta algorithm

For the implementation of the explicit Runge-Kutta steppers a new method using template metaprogramming is used. With this method a Runge-Kutta algorithm is only specified by its coefficients. Of course, it is also possible to use classical programming techniques for this task, but one will certainly use performance. The template metaprogramming algorithm avoids this downfall. Only compile-time container and algorithms are used here, such that the compiler can perform powerful optimization. The performance of these algorithms is comparable to its hand-written version, see also the next section. The complete algorithm is described in detail here (in German only) [5], and here [7].

Performance

The performance of odeint has been carefully checked against existing solutions like the methods of the GNU Scientific library (gsl) or Numerical Recipes (NR). It turned out that the abstraction introduced in odeint does not affect its performance assuming that a modern C++ compiler with optimization is used. Performance result for the Lorenz system and boost::array as state type are shown in the figure below. The percentage values are related to the first version of odeint. odeint stands for the first version of odeint, generic for the current implementation in odeint v2, NR for the methods from the Numerical Recipes, gsl for the GNU Scientific library and rt_gen for a generic run-time implementation of the method under consideration.

Summary and outlook

In this article odeint v2 has been introduced - a high-level C++ library for solving ordinary differential equations. Its main enhancement over the first version are the use of algebras and operations, two concepts which allow one to change the way the basic operations in odeint are performed. Using these techniques one can solve ODEs on modern graphic cards, or one can use compile-time containers and advanced dimension-analysis classes like Boost.Fusion and Boost.Units. Furthermore, more numerical methods have been implemented and dense-output functionality has been introduced.

The following table shows all steppers which are currently implemented in odeint

Method

Class name

Comment

Explicit Euler

euler

Very simple, only for demonstrating purpose

Runge-Kutta 4

runge_kutta4

The classical Runge Kutta scheme, good general scheme without error control

Cash-Karp

runge_kutta_cash_karp54

Good general scheme with error estimation

Dormand-Prince 5

runge_kutta_dopri5

Standard method with error control and dense output

Fehlberg 78

runge_kutta_fehlberg78

Good high order method with error estimation

Adams-Bashforth-Moulton

adams_bashforth_moulton

Multi-step method with high performance

Controlled Error Stepper

controlled_runge_kutta

Error control for the Runge-Kutta steppers

Dense Output Stepper

dense_output_runge_kutta

Dense output for the Runge-Kutta steppers

Bulirsch-Stoer

bulirsch_stoer

Stepper with step size, order control and dense output. Very good if high precision is required.

Implicit Euler

implicit_euler

Basic implicit routine

Rosenbrock 4

rosenbrock4

Solver for stiff systems with error control and dense output

Symplectic Euler

symplectic_euler

Basic symplectic solver for separable Hamiltonian system

Symplectic RKN McLachlan

symplectic_rkn_sb3a_mclachlan

Symplectic solver for separable Hamiltonian system with order 6

odeint is currently under heavy development. Therefore, some class names and include files might still change and new features and methods might be implemented. The source code provided with this article is a snapshot of the current development.

If you like to contribute to this library, or have some interesting points, criticism, or suggestions, you are cordially invited.

About the Author

Comments and Discussions

I am simulating an electrical circuit using odeint. The results are not compatible with matlab results (using ode23). I have tried fixed step size adaptive stpe size and dense output. I have tried reducing the step sizes, and error tolerances. I know that matlab results are correct because I simulated my circuit in a dedicated circuit simulator, and matlab matches that.

Would you have any suggestions on debugging the problem? If needed, I can provide my code, which is not large.

Many many thanks. Exactly what I needed. I'm not very expert about c++ so I didn't understand how really Thrust works with Odeint. Now I did a step forward.
How about including an example like this in the odeint tutorial? It's already a very good work but given that all the examples are about autonomous systems, it can be more effective with a simple example like this.
Thanks,
Michele

Hi,
sorry to disturb but I am new to the code project and excited about using odeint v2 for academic research. I am trying to use the implicit euler stepper for the first step in a system of nonlinear ODEs and when I type stepper.do_step (system, x, r , dr) I get an error message saying the the function do_step can not be found. I thought that this was a function every stepper had. Am I doing something wrong? I tried looking for an example with the implicit euler used but I have not seen any.
Any help would be appreciated. Thank you in advance!

maybe you are misusing the stepper method of the implict_euler. As system function it expects a pair of function object, one returning the r.h.s. of the ode and the other one the Jacobian of the r.h.s. An example can be found here (ok, it is the rosenbrock stepper, but this one is very similar to the implicit euler).

Hi,
thank you for your quick reply. I finally got to try the example you gave me with implicit euler but I still have problems. To be honest, I'm pretty new to template usage and probably my mistakes are trivial but I still need help. I just need to do the first step in my problem implicitly and use runge kutta for the rest. Below I have written the code I have so far, using the stiff system from your web page. Please, let me know what I am doing wrong!

1. the implicit_euler is instantiated only with one template argument double
2. the Jacobian is calculated without explicitly calculating the partial derivative of the system function with respect to the time. Simply replace stiff_system_jacobi with

There are some small inconsistencies with the implicit steppers. The Rosenbrock stepper (another implicit stepper) expect the partial derivative of the system function w.r.t. the time explicitly, while the implicit Euler does not.

Hi,
thank you for your help. I opted to use the Rosenbrock method since my system of ODE has explicit t dependence. However, I still have problems. This time it seems to be due to division by t^2 at t=0 which is the beginning of my domain. I thought that since it is an implicit method it should deal wtih this kind of situations. The message I get is "This application has requested the Runtime to terminate it in an unusual way. Please Contact the application's support team for more information." And my code is below:

Hi,
I think I already mentioned that I'm starting at 0. But that should be a problem only for explicit methods whcih use the current point to go forward. Impicit methods like backward Euler should use the t+dt point to find the value of the function at t+dt, which is the main reason why I want to use your code. I thought Rosenbrock had the same capbility. Am I wrong? Is there a way to do that? I have already tried using 0.01 instead of 0 but that leads to big numerical errors later in calculations. Thank you!

Ok, Rosenbrock is a semi-implicit method. It calculates the Jacobian at t=0 and evaluates the function at t=0. I think you can use implicit euler method. But it expects the system function in a slightly different way, such that you have to rewrite it. But this should only a minor change for your problem.

Hi,
I did the changes you have previously mentioned but neither do step nor integrate_adaptive seem to work with the Euler method. You mentioned that I have to rewrite the function. Would you mind showing me how? There's no example on how to use implicit_euler. My changed code is below:

Note, that the implicit_euler does not have any dense output or step size control, such that integrate_adaptive is only an integrate_const. Sorry for the inconsistency of the system functions interface. Unifying the implicit steppers system function is an important point on our ToDo list. Of course, we also plan to implement more implicit methods, but it will take some time.

Hi,
I used the implicit euler stepper together with runge kutta 4 and they worked like a charm! Thank you for the help provided. By the way, I think it would be a good idea to add examples on how each stepper works since there some differences in their implementations. For instance, integrate-adaptive shouldn't be used at all with implicit Euler since there is integrate_const which gives better results. I might bug you again in the near future!

Ok, thank you for the hint. We will add an example about the implicit_euler in the tutorial section.

By the way, in the case of the implicit_euler integrate_adaptive and integrate_const should give the same results. Of course, integrate_adaptive is misleading since it simply redirects in this case to integrate_adaptive.

yesterday i downloaded the odeint-v2 package, but i have no clue how to compile or setup it.
I have a pc with linux (opensuse 12.1) and boost 1.46-1.

I know the classic unix-tools like make, gcc and friends, (also c and some c++) but i do not know how to use bjam,
bootstrap and this. i googled around, reading the boost site but this did not help me.

Can one friendly soul give me a hint how to start to setup odeint-v2 from zero?

There is nothing to setup. odeint is completely header only, all you have to do is to specify the directory during compilation. For gcc this is easy, just add the compiler flags -I/path/to/odeint -I/path/to_boost.

If you want to compile the examples and unit tests of odeint you have to

export BOOST_ROOT=/usr/share/doc/packages/boost-1.46.1 ((Here boost is installed when using the opensuse dvd))

# ~/downloads/odeint-v2> bjam
Unable to load Boost.Build: could not find "boost-build.jam"
---------------------------------------------------------------
Attempted search from /home/mma/downloads/odeint-v2 up to the root
at /usr/share/boost-build
and in these directories from BOOST_BUILD_PATH and BOOST_ROOT: /usr/share/boost-build, /usr/share/doc/packages/boost-1.46.1.
Please consult the documentation at 'http://www.boost.org'.

i fear there is lacking something in the installation. find only finds:

i succeeded to setup odeint-v2 and could start all example
But when i tried to compile the thrust examples for using cuda it failed.
(I have installed cuda 4.1, gcc3.6, thrust 1.48.0 and 2 nvidia gtx580 cards):

well because i am interested in solving odes on a gpu i was interested to see how this works with "thrust"
as i understand from your article thrust pushed and pulls the code (ordinarily spoken) from and to the gpu
so i wanted to see how this works.

watching the gpu indead show that she is idle most of the time.
but the cpu (i7-2600 4 cores with 2 threads) show that all 8 virtual cpus are working fullspeed)

So, you have a speedup factor of approx. 10. This looks ok. In the meanwhile we found some bug in thrust_resizer which might be responsible for the crashes of dopri5 in the phase_oscillator_ensemble example. I am not sure if it is fixed right now, but will keep you informed.

i just commited a bugfix for the thrust backend in odeint. i hope this helps with your segmentation faults. sorry for this, i forgot to adjust the thrust backend when applying some changes in the memory management.

Sorry for the late reply. Solving a very simple ODE does not make sense with CUDA, unless you are solve many of them in parallel. What kind of ODE are you familar with? I could provide an example for dxdt = -a x, where several ODEs with different parameters a are solved in parallel. Is this ok for you?