How to perform simple numeric integration with odeint in C++

Question

Can you provide me with a simple example of performing a numeric integration with odeint in C++?

I would like to use the convenient integrate function, documented as:

integrate( system , x0 , t0 , t1 , dt )

Also I'm not sure, how to pass it instead of a function or a functor, a class method, if that's possible.

Answer 1

In C++11 you can use a simple lambda function wrapping the call to your member method

Class c;
auto f = [&c]( const state_type & x , state_type &dxdt , double t ) {
    c.system_func( x , dxdt , t ); };
integrate( f , x0 , t0 , t1 , dt );

std::bind might also work, but then you have to take care if values are passed by reference of by value.

In C++03 you need to write a simple wrapper around your class method

struct wrapper
{
    Class &c;
    wrapper( Class &c_ ) : c( c_ )  { }
    template< typename State , typename Time >
    void operator()( State const &x , State &dxdt , Time t ) const
    {
        c.system_func( x , dxdt , t );
    }
};

// ...
integrate( wrapper( c ) , x0 , t0 , t1 , dt );

(Boost.Bind will not work correctly with more then two arguments).

Answer 2

You mean examples in addition to the ones provided online?

#include <iostream>
#include <boost/array.hpp>

#include <boost/numeric/odeint.hpp>

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

const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;

typedef boost::array< double , 3 > state_type;

void lorenz( const state_type &x , state_type &dxdt , double t )
{
    dxdt[0] = sigma * ( x[1] - x[0] );
    dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
    dxdt[2] = -b * x[2] + x[0] * x[1];
}

void write_lorenz( const state_type &x , const double t )
{
    cout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] << endl;
}

int main(int argc, char **argv)
{
    state_type x = { 10.0 , 1.0 , 1.0 }; // initial conditions
    integrate( lorenz , x , 0.0 , 25.0 , 0.1 , write_lorenz );
}

And regarding the system, you can provide anything where the following expression is valid:

sys( x , dxdt , t ) // returning void

Check the user's guide (and more examples) online.