Boost C++ Libraries

As a first example we take the most simple case and implement our own vector my_vector which will provide a container interface. This makes Boost.Range working out-of-box. We add a little functionality to our vector which makes it allocate some default capacity by construction. This is helpful when using resizing as then a resize can be assured to not require a new allocation.

template< size_t MAX_N >
class my_vector
{
    typedef std::vector< double > vector;

public:
    typedef vector::iterator iterator;
    typedef vector::const_iterator const_iterator;

public:
    my_vector( const size_t N )
        : m_v( N )
    {
        m_v.reserve( MAX_N );
    }

    my_vector()
        : m_v()
    {
        m_v.reserve( MAX_N );
    }

// ... [ implement container interface ]

The only thing that has to be done other than defining is thus declaring my_vector as resizeable:

// define my_vector as resizeable

namespace boost { namespace numeric { namespace odeint {

template<size_t N>
struct is_resizeable< my_vector<N> >
{
    typedef boost::true_type type;
    static const bool value = type::value;
};

} } }

If we wouldn't specialize the is_resizeable template, the code would still compile but odeint would not adjust the size of temporary internal instances of my_vector and hence try to fill zero-sized vectors resulting in segmentation faults! The full example can be found in my_vector.cpp

If your state type does work with Boost.Range, but handles resizing differently you are required to specialize two implementations used by odeint to check a state's size and to resize:

As an example we will use a std::list as state type in odeint. Because std::list is not supported by boost::size we have to replace the same_size and resize implementation to get list to work with odeint. The following code shows the required template specializations:

typedef std::list< double > state_type;

namespace boost { namespace numeric { namespace odeint {

template< >
struct is_resizeable< state_type >
{ // declare resizeability
    typedef boost::true_type type;
    const static bool value = type::value;
};

template< >
struct same_size_impl< state_type , state_type >
{ // define how to check size
    static bool same_size( const state_type &v1 ,
                           const state_type &v2 )
    {
        return v1.size() == v2.size();
    }
};

template< >
struct resize_impl< state_type , state_type >
{ // define how to resize
    static void resize( state_type &v1 ,
                        const state_type &v2 )
    {
        v1.resize( v2.size() );
    }
};

} } }

With these definitions odeint knows how to resize std::lists and so they can be used as state types. A complete example can be found in list_lattice.cpp.

In the following example we will try to use the gsl_vector type from GSL (GNU Scientific Library) as state type in odeint. We will realize this by implementing a wrapper around the gsl_vector that takes care of construction/destruction. Also, Boost.Range is extended such that it works with gsl_vectors as well which required also the implementation of a new gsl_iterator.

	Note
	odeint already includes all the code presented here, see gsl_wrapper.hpp, so gsl_vectors can be used straight out-of-box. The following description is just for educational purpose.

The GSL is a C library, so gsl_vector has neither constructor, nor destructor or any begin or end function, no iterators at all. So to make it work with odeint plenty of things have to be implemented. Note that all of the work shown here is already included in odeint, so using gsl_vectors in odeint doesn't require any further adjustments. We present it here just as an educational example. We start with defining appropriate constructors and destructors. This is done by specializing the state_wrapper for gsl_vector. State wrappers are used by the steppers internally to create and manage temporary instances of state types:

template<>
struct state_wrapper< gsl_vector* >
{
    typedef double value_type;
    typedef gsl_vector* state_type;
    typedef state_wrapper< gsl_vector* > state_wrapper_type;

    state_type m_v;

    state_wrapper( )
    {
        m_v = gsl_vector_alloc( 1 );
    }

    state_wrapper( const state_wrapper_type &x )
    {
        resize( m_v , x.m_v );
        gsl_vector_memcpy( m_v , x.m_v );
    }


    ~state_wrapper()
    {
        gsl_vector_free( m_v );
    }

};

This state_wrapper specialization tells odeint how gsl_vectors are created, copied and destroyed. Next we need resizing, this is required because gsl_vectors are dynamically sized objects:

template<>
struct is_resizeable< gsl_vector* >
{
    typedef boost::true_type type;
    const static bool value = type::value;
};

template <>
struct same_size_impl< gsl_vector* , gsl_vector* >
{
    static bool same_size( const gsl_vector* x , const gsl_vector* y )
    {
        return x->size == y->size;
    }
};

template <>
struct resize_impl< gsl_vector* , gsl_vector* >
{
    static void resize( gsl_vector* x , const gsl_vector* y )
    {
        gsl_vector_free( x );
        x = gsl_vector_alloc( y->size );
    }
};

Up to now, we defined creation/destruction and resizing, but gsl_vectors also don't support iterators, so we first implement a gsl iterator:

/*
 * defines an iterator for gsl_vector
 */
class gsl_vector_iterator
      : public boost::iterator_facade< gsl_vector_iterator , double ,
                                       boost::random_access_traversal_tag >
{
public :

    gsl_vector_iterator( void ): m_p(0) , m_stride( 0 ) { }
    explicit gsl_vector_iterator( gsl_vector *p ) : m_p( p->data ) , m_stride( p->stride ) { }
    friend gsl_vector_iterator end_iterator( gsl_vector * );

private :

    friend class boost::iterator_core_access;
    friend class const_gsl_vector_iterator;

    void increment( void ) { m_p += m_stride; }
    void decrement( void ) { m_p -= m_stride; }
    void advance( ptrdiff_t n ) { m_p += n*m_stride; }
    bool equal( const gsl_vector_iterator &other ) const { return this->m_p == other.m_p; }
    bool equal( const const_gsl_vector_iterator &other ) const;
    double& dereference( void ) const { return *m_p; }

    double *m_p;
    size_t m_stride;
};

A similar class exists for the const version of the iterator. Then we have a function returning the end iterator (similarly for const again):

gsl_vector_iterator end_iterator( gsl_vector *x )
{
    gsl_vector_iterator iter( x );
    iter.m_p += iter.m_stride * x->size;
    return iter;
}

Finally, the bindings for Boost.Range are added:

// template<>
inline gsl_vector_iterator range_begin( gsl_vector *x )
{
    return gsl_vector_iterator( x );
}

// template<>
inline gsl_vector_iterator range_end( gsl_vector *x )
{
    return end_iterator( x );
}

Again with similar definitions for the const versions. This eventually makes odeint work with gsl vectors as state types. The full code for these bindings is found in gsl_wrapper.hpp. It might look rather complicated but keep in mind that gsl is a pre-compiled C library.

As seen above, the standard way of performing algebraic operations on container-like state types in odeint is to iterate through the elements of the container and perform the operations element-wise on the underlying value type. This is realized by means of the range_algebra that uses Boost.Range for obtaining iterators of the state types. However, there are other ways to implement the algebraic operations on containers, one of which is defining the addition/multiplication operators for the containers directly and then using the vector_space_algebra. If you use this algebra, the following operators have to be defined for the state_type:

Defining these operators makes your state type work with any basic Runge-Kutta stepper. However, if you want to use step-size control, some more functionality is required. Specifically, operations like max_i( |err_i| / (alpha * |s_i|) ) have to be performed. err and s are state_types, alpha is a scalar. As you can see, we need element wise absolute value and division as well as an reduce operation to get the maximum value. So for controlled steppers the following things have to be implemented:

Here we show how to implement the required operators on a state type. As example we define a new class point3D representing a three-dimensional vector with components x,y,z and define addition and scalar multiplication operators for it. We use Boost.Operators to reduce the amount of code to be written. The class for the point type looks as follows:

class point3D :
    boost::additive1< point3D ,
    boost::additive2< point3D , double ,
    boost::multiplicative2< point3D , double > > >
{
public:

    double x , y , z;

    point3D()
        : x( 0.0 ) , y( 0.0 ) , z( 0.0 )
    { }

    point3D( const double val )
        : x( val ) , y( val ) , z( val )
    { }

    point3D( const double _x , const double _y , const double _z  )
        : x( _x ) , y( _y ) , z( _z )
    { }

    point3D& operator+=( const point3D &p )
    {
        x += p.x; y += p.y; z += p.z;
        return *this;
    }

    point3D& operator*=( const double a )
    {
        x *= a; y *= a; z *= a;
        return *this;
    }

};

By deriving from Boost.Operators classes we don't have to define outer class operators like operator+( point3D , point3D ) because that is taken care of by the operators library. Note that for simple Runge-Kutta schemes (like runge_kutta4) only the + and * operators are required. If, however, a controlled stepper is used one also needs to specify the division operator / because calculation of the error term involves an element wise division of the state types. Additionally, controlled steppers require an abs function calculating the element-wise absolute value for the state type:

// only required for steppers with error control
point3D operator/( const point3D &p1 , const point3D &p2 )
{
    return point3D( p1.x/p2.x , p1.y/p2.y , p1.z/p2.z );
}

point3D abs( const point3D &p )
{
    return point3D( std::abs(p.x) , std::abs(p.y) , std::abs(p.z) );
}

Finally, we have to provide a specialization to calculate the infintity norm of a state:

// also only for steppers with error control
namespace boost { namespace numeric { namespace odeint {
template<>
struct vector_space_norm_inf< point3D >
{
    typedef double result_type;
    double operator()( const point3D &p ) const
    {
        using std::max;
        using std::abs;
        return max( max( abs( p.x ) , abs( p.y ) ) , abs( p.z ) );
    }
};
} } }

Again, note that the two last steps were only required if you want to use controlled steppers. For simple steppers definition of the simple += and *= operators are sufficient. Having defined such a point type, we can easily perform the integration on a Lorenz system by explicitely configuring the vector_space_algebra in the stepper's template argument list:

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

void lorenz( const point3D &x , point3D &dxdt , const double t )
{
    dxdt.x = sigma * ( x.y - x.x );
    dxdt.y = R * x.x - x.y - x.x * x.z;
    dxdt.z = -b * x.z + x.x * x.y;
}

using namespace boost::numeric::odeint;

int main()
{

    point3D x( 10.0 , 5.0 , 5.0 );
    // point type defines it's own operators -> use vector_space_algebra !
    typedef runge_kutta_dopri5< point3D , double , point3D ,
                                double , vector_space_algebra > stepper;
    int steps = integrate_adaptive( make_controlled<stepper>( 1E-10 , 1E-10 ) , lorenz , x ,
                                    0.0 , 10.0 , 0.1 );
    std::cout << x << std::endl;
    std::cout << "steps: " << steps << std::endl;
}

The whole example can be found in lorenz_point.cpp

	Note
	For the most state_types, odeint is able to automatically determine the correct algebra and operations. But if you want to use your own state_type, as in this example with point3D, you have to manually configure the right algebra/operations, unless your state_type works with the default choice of range_algebra and default_operations.