## Numerical integration in C++ with Boost odeint

May 2016, Bruxelles

Back to HOMEPAGE

In this notes there's nothing you can't find in the official documentation, moreover informations here may be misleading or even wrong. I've written this page just to save some time to whoever needs a quick kickstart to numerical integration with the Boost library (included myself). I'll appreciate any comment, suggestion and critics!

A very common task in numerical computing is integrating systems of differential equations. The boost library makes it very easy, but understanding how to do it the first time by searching into the official documentation and surfing the examples is time consuming, so I've written this brief howto.
This howto was written under GNU/Linux with Boost v1.60.

Some references: And here is the index:

INTRO:
Setting up the Boost C++ library
TUTORIAL:
0 - Integration functions
1 - Simplest ODE problem
2 - Integration functions: overview
3 - Integration with steppers
4 - Stiff problems and Jacobian

#### Setting up the Boost C++ library

Good news: the Boost library implements many functions into header files, so in many cases no compilation is necessary! This is our case. Download the boost library (version higher than 1.53 since odeint was introduced in that Boost version) from www.boost.org, and extract it where you like.
In my case, I've extracted the Boost library into the directory: ``` ~/boost_1_60_0/ ```
After extracting, you should see the directory boost/numeric, in which there's the odeint.hpp file. We will use this file including it into our C++ source.

#### 0 - Integration functions

In case you are familiar to the Matlab/Octave syntax for the ode45 integrators, using odeint is quite similar.
First of all, odeint have 5 integration functions:
1. integrate_const
2. integrate_n_steps
3. integrate_times
5. integrate
The syntax for integrate_const is (the others vary just a little, see below):
` integrate_const( stepper, system, x0, t0, t1, dt, observer)`

Where:
• stepper is the mathematical scheme you wish to use (such as Runge-Kutta or Euler). It defines how the integration is to be made: constand steps, adaptive, stiff systems etc
• system is a function implementing the mathematical system (or the scalar equation) you want to solve. This function takes as an input the current state and time and creates the time derivative dxdt
• x0 is the initial condition
• t0 and t1 are initial and final times for integration
• dt is a timestep, but its signifiance depends on the type of stepper
• observer this is a function that is called every N timesteps, that you can use to print intermediate solutions. It's optional. Without it the integrator will just integrate from t0 to t1 and then return with no useful output

#### 1 - Simplest ODE problem

The following code defines a scalar ordinary differential equation and solves it using the simplest integration function.
``` #include <iostream>
#include <boost/numeric/odeint.hpp>       // odeint function definitions

using namespace std;
using namespace boost::numeric::odeint;

// Defining a shorthand for the type of the mathematical state
typedef std::vector< double > state_type;

// System to be solved: dx/dt = -2 x
void my_system( const state_type &x , state_type &dxdt , const double t )
{    dxdt = -2*x;   }

// Observer, prints time and state when called (during integration)
void my_observer( const state_type &x, const double t )
{    std::cout  << t << "   " << x << std::endl;   }

// ------  Main
int main()
{
state_type x0(1); // Initial condition, vector of 1 element (scalar problem)
x0 = 10.0;

// Integration parameters
double t0 = 0.0;
double t1 = 10.0;
double dt = 0.1;

// Run integrator
integrate( my_system, x0, t0, t1, dt, my_observer );
}
```
In order to compile you have to include the path for the extracted Boost folder with the -I option. For example, my Boost directory is ~/boost_1_60_0/ so I run: BOOSTLIB=~/boost_1_60_0 g++ -I\$BOOSTLIB my_file.cpp The output of the program should be something like: 0 10 0.1 8.18731 0.222349 6.41018 0.344698 5.0188 0.467046 3.92942 0.589395 3.07651 0.711744 2.40873 0.838353 1.86989 0.964962 1.45159 1.09157 1.12687 1.2236 0.865356 [...] A couple of notes about the program.
1. First of all, as you can see the integrate function does not require a stepper. In fact, integrate is the simplest integration function, used just to perform some integration on the fly. You usually want to use the others (see below).
2. Secondly, you should note that the solution of the mathematical problem, x (aka the "state" of the system) have a type that we called state_type. One can choose various different types, here I've opted for a simple std vector of doubles, whereas we could have used ublas types for example.
3. The my_system function takes as an input the current state and fills the second argument, dxdt that is the left hand side of the ODE problem.
4. The observer just prints the intermediate result. We could have used it to write the results on a file.
5. The observer is invoked roughly every "dt".
Comparison with exact solution: #### 2 - Integration functions: overview

The choice of the integration function to be used mainly depends on what kind of output you need, that is how often you need to call the observer.
• integrate_const calls the observer at equidistant times separated by dt. The syntax is:
` integrate_const( stepper , system , x0 , t0 , t1 , dt , observer ) `
• integrate_n_steps is like the previous one but does not ask for the end time, instead it asks for the number of timesteps to integrate. Syntax:
` integrate_n_steps( stepper , system , x0 , t0 , dt , n , observer ) `
• integrate_times performs integration at given timesteps. Syntax:
` integrate_times( stepper , system , x0 , times_start , times_end , dt , observer )` ` integrate_times( stepper , system , x0 , time_range , dt , observer ) `
• integrate_adaptive is used if you want to call the observer at each timestep. Syntax:
` integrate_adaptive( stepper , system , x0 , t0 , t1 , dt , observer ) `
• integrate this is the "convenience integrate function", simple and quick to set without need to specify a stepper. Syntax:
` integrate( system , x0 , t0 , t1 , dt , observer ) `

Note that the choice of the stepper (think to an adaptive method, that changes the integration timestep according to the error) will modify how actually the integration is performed, but calls to the observer (and thus your output) follows the rules above.
Steppers are the actual "solvers", and are called recursively by the choosen integration function.

#### 3 - Integration with steppers

Let's now switch to the integrate_const function and modify the previous program to introduce a stepper. With a stepper we can choose the type of integration to be performed, such as Runge-Kutta method, adaptive method or whatever.
For example, to choose a 4th order Runge-Kutta method you can call the integration function like this:
``` integrate_const( runge_kutta4<state_type>(), my_system, x0, t0, t1, dt, my_observer );
```

Or you might prefere a cleaner syntax:
``` typedef runge_kutta4<state_type> rk4;
integrate_const( rk4(), my_system, x0, t0, t1, dt, my_observer );
```

Steppers can also be used to directly perform an integration step, in the following way. Here only the main() is shown, the rest being the very same as above:
``` // ------  Main
int main()
{
state_type x(1); // the state, vector of 1 element (scalar problem)
x = 10.0;

// Integration parameters
double t = 0.0;
double dt = 0.1;
size_t nSteps = 1000;

// Create stepper
runge_kutta4<state_type> rk4;

// Perform integration step by step
for(int ii = 0; ii < nSteps; ++ii)
{
t += dt;

// now perform only one integration step. Solution x is overwritten
rk4.do_step(my_system, x, t, dt);

// ... and output results
cout << t << "   " << x << endl;
}
}
```
Note that the x variable is initialized at the initial condition and is rewritten and output at each step.

As already stated, steppers are basically "solvers". There are some different kinds of steppers:
• basic stepper
• error stepper
• controlled stepper
• dense output stepper
The basic stepper have already been used and other possible integration methods are: runge_kutta4, euler, runge_kutta_cash_karp54, runge_kutta_dopri5, runge_kutta_fehlberg78, modified_midpoint, rosenbrock4.
Error steppers are steppers that provide an error estimation. Those steppers are required to build controlled steppers, that can adapt the integration timestep. Some error steppers are runge_kutta_cash_karp54, runge_kutta_dopri5 and runge_kutta_fehlberg78.
Controlled steppers are built on error steppers. With controlled steppers, the routine starts performing the integration with the given timestep size dt, but may decide with an error criteria to use a different timestep in case the suggested dt was not adequate. To build a controlled stepper for Runge Kutta, you can use controlled_runge_kutta.
Dense steppers perform some kind of continuous interpolation between timesteps. A dense stepper is the rosenbrock4, useful for stiff systems (see below).

Example
In the following example we'll use the controlled stepper runge_kutta_cash_karp54 to solve a slightly more interesting problem: It's the equation of the harmonic oscillator. For sure you know that with a choice of ξ > 1 the equation shows a nonoscillating behaviour (supercritical damping). By choosing the appropriate initial conditions, the solution can present an overshoot. Here is the equation that we are going to solve plus BCs, along with a sketch of the exact solution. Note that since odeint deals with first order equations, the second order harmonic oscillator have to be recasted into a second order ODE system: In the following code I define a dimension-2 linear system and solve it using two different steppers. The first one is the already seen runge_kutta4 (basic stepper), the latter being a controller stepper based on the runge_kutta_cash_karp54 (that internally adjusts the integration steps based on error considerations).
The timestep is set very big on purpose, so that the runge_kutta4 basic stepper fails, while the adaptive stepper obtains a good result, internally automatically refining the timestep.

Take a look at the code and don't miss the notes right after!
```#include <iostream>

#include <boost/numeric/odeint.hpp> // odeint function definitions

using namespace std;
using namespace boost::numeric::odeint;

// Defining a shorthand for the type of the mathematical state
typedef std::vector< double > state_type;

// System to be solved: dx/dt = -2 x
void my_system( const state_type &x , state_type &dxdt , const double t )
{
dxdt =  0*x  +   1*x;
dxdt = -1*x  - 2.2*x;
}

// Observer, prints time and state when called (during integration)
void my_observer( const state_type &x, const double t )
{
std::cout  << t << "   " << x << "   " << x << std::endl;
}

// ------  Main
int main()
{
state_type x0(2); // Initial condition, vector of 2 elements (position and velocity)
x0 = 0.0;
x0 = 1.0;

// Integration parameters
double t0 = 0.0;
double t1 = 20.0;
double dt = 1.0;

// ----  Steppers definition  ----
// Basic stepper:
// follows given timestep size "dt"
typedef runge_kutta4<state_type> rk4;

// Error stepper, used to create the controlled stepper
typedef runge_kutta_cash_karp54< state_type > rkck54;

// Controlled stepper:
// it's built on an error stepper and allows us to have the output at each
// internally defined (refined) timestep, via integrate_adaptive call
typedef controlled_runge_kutta< rkck54 > ctrl_rkck54;

// ----  Run integrations with the different steppers  ----

// Run integrator with rk4 stepper
std::cout << "==========  rk4 - basic stepper  ====================" << std::endl;
integrate_adaptive( rk4(), my_system, x0, t0, t1, dt, my_observer );

// Run integrator with controlled stepper
std::cout << "==========  ctrl_rkck54 - Controlled Stepper  =======" << std::endl;
x0 = 0.0; // Reset initial conditions
x0 = 1.0;
integrate_adaptive( ctrl_rkck54(), my_system, x0 , t0 , t1 , dt, my_observer );
}
```

In solid line the exact solution, circles and squares are numerical solutions. A couple of notes:
• I've choosed the integrate_adaptive integration function, that calls the solution also at intermediate timesteps, so that we can see how the controlled stepper works.
• When a controlled stepper is passed to integrate_adaptive, the timestep dt is just an initial guess timestep. This timestep as well as all the others can be refined with error criteria, in fact you can see from the plot at right that despite the timestep was initially set to dt = 1.0, at the end of the controlled stepper integration the actual step is bigger!
Here is the output of the program:
```==========  rk4 - basic stepper  ====================
0   0   1
1   0.279667   0.0914
2   0.223193   -0.0698595
3   0.138185   -0.0688047
4   0.0784087   -0.0449346
5   0.0428421   -0.0260353
6   0.0229939   -0.0143611
7   0.0122327   -0.00774322
8   0.00647889   -0.0041288
9   0.00342373   -0.0021893
10   0.00180716   -0.00115761
11   0.000953317   -0.000611208
12   0.000502743   -0.000322475
13   0.000265086   -0.000170075
14   0.000139763   -8.96806e-05
15   7.36854e-05   -4.72839e-05
16   3.88473e-05   -2.49291e-05
17   2.04802e-05   -1.31428e-05
18   1.07971e-05   -6.92889e-06
19   5.69216e-06   -3.65289e-06
20   3.00087e-06   -1.92578e-06
==========  ctrl_rkck54 - Controlled Stepper  =======
0   0   1
0.2   0.160729   0.62909
0.4   0.25906   0.369921
0.6   0.314033   0.191075
0.809566   0.339954   0.065062
1.01913   0.344383   -0.0166857
1.2428   0.33413   -0.0702323
1.47763   0.31359   -0.101238
1.73025   0.285836   -0.11597
2.00191   0.253734   -0.118652
2.29508   0.219629   -0.112967
2.6131   0.185377   -0.101919
2.96069   0.152372   -0.087867
3.34299   0.121725   -0.0726501
3.76712   0.0941826   -0.0576185
4.24293   0.0702024   -0.0437071
4.78465   0.049993   -0.0315045
5.41432   0.0335592   -0.0213196
6.17034   0.0207314   -0.0132375
7.13966   0.0111524   -0.00714214
8.2138   0.0056026   -0.0035926
9.33462   0.00272999   -0.00175145
10.5825   0.0012258   -0.00078657
11.998   0.00049426   -0.000317177
13.6366   0.000172728   -0.000110845
15.5745   4.98493e-05   -3.199e-05
17.9281   1.10743e-05   -7.10652e-06
20   2.9352e-06   -1.8835e-06
```

Generation functions
Odeint provides two handy functions to generate controlled and dense_output steppers, called make_controlled and make_dense_output. Here is the documentation page for Boost version 1.60.
Here's an example:
```// Define the error stepper
typedef runge_kutta_cash_karp54< state_type > error_stepper_rkck54;

// Error bounds
double err_abs = 1.0e-10;
double err_rel = 1.0e-6;

// Fire!
integrate_adaptive( make_controlled( err_abs , err_rel , error_stepper_rkck54() ) ,
my_system , x0 , t0 , t1 , dt , my_observer );
```

#### 4 - Stiff problems and Jacobian

Stiff problems are characterized by very different timescales. Think for example at a chemically reacting mixture: chemical rates may be very different among each others (by orders of magnitude!) That characteristic makes common integrators inappropriate and one needs apposite solvers.
Consider for example the following system: This is basically a system where the derivative of the first component x1 is likely way bigger than the one for the second component x2. A quick check with matlab shows that a Runge-Kutta 45 adaptive method (ode45) takes 70000 steps to solve it into a certain interval, whereas a stiff solver (ode15) takes only 100 timesteps!
The more the system stiffness, the more the disparity: what if we increase the stiffness, then? Let's put the coefficients of the first equation to 30000 and 20000, the required timesteps for the nonstiff solver increase to 700000, while the stiff solver takes still 100 timesteps!!! Of course, nobody knows what kind of black magic they put into Matlab functions, but I think you are now convinced of the importance of stiff solvers.

Odeint implements some stiff solvers, such as the rosenbrock4, but they additionaly require to provide the system Jacobian, that is the matrix of the partial derivatives. Also, note that the stiff integrator of Odeint *requires* that the system state is of ublas::vector type and the jacobian *must* be a ublas::matrix!

Here's a program that solves the problem defined above. Note that:

• I'm using a slightly different syntax for the system this time, not defining it as a function, but using functors: the system is defined as a structure and the operator "()" is overloaded so that when it is called with parenthesis it behaves like a function.
• The Jacobian too is defined as a struct, in the very same way as my_system
• I've used integrate_const, that differs from the previously used integrate_adaptive only in that we obtain the output precisely at each "dt"
• In the observer I use the std function setw(), that gives a slightly more cute output.
``` #include <iostream>

#include <boost/numeric/odeint.hpp> // odeint function definitions

using namespace std;
using namespace boost::numeric::odeint;

// Defining shorthands for boost types vector and matrix.
// vector is used for the system state "x" and matrix for the jacobian "J"
typedef boost::numeric::ublas::vector< double > vector_type;
typedef boost::numeric::ublas::matrix< double > matrix_type;

// System to be solved.
// This is a functor: "operator()" overloads the operator "()" making the
// struct behave like a function
struct my_system
{
void operator()( const vector_type &x, vector_type &dxdt, const double t)
{
dxdt = -3000*x -2000*x;
dxdt =    -2*x    -3*x;
}
};

// Jacobian of the system.
// Functor! See struct "my_system"
struct jacobian
{
void operator()( const vector_type &x, matrix_type &J, double t, vector_type &dfdt )
{
J(0,0) = -3000;
J(0,1) = -2000;
J(1,0) = -2;
J(1,1) = -3;
}
};

// Observer, prints time and state when called (during integration)
void my_observer( const vector_type &x, const double t )
{
std::cout << t << setw(15) << x << setw(15) << x << std::endl;
}

// ------  Main
int main()
{
vector_type x0(2); // Initial condition, vector of 2 elements (position and velocity)
x0 = 0.0;
x0 = 1.0;

// Integration parameters
double t0 = 0.0;
double t1 = 20.0;
double dt = 1.0;

// ----  Run integrations with the different steppers  ----

// Stiff integration
integrate_const(make_dense_output< rosenbrock4< double > >(1.0e-6, 1.0e-6),
make_pair(my_system(), jacobian()),
x0, t0, t1, dt, my_observer);
}
```

That's it. I hope you enjoyed!
FOR ANY CRITICS/SUGGESTIONS/WHATEVER CONTACT ME!
For any spotted typo I will send you a cookie.

Cheers!

Back to HOMEPAGE