...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
Home  Libraries  People  FAQ  More 
In odeint the stepper algorithms are implemented independently of the underlying
fundamental mathematical operations. This is realized by giving the user
full control over the state type and the mathematical operations for this
state type. Technically, this is done by introducing three concepts: StateType,
Algebra, Operations. Most of the steppers in odeint expect three class types
fulfilling these concepts as template parameters. Note that these concepts
are not fully independent of each other but rather a valid combination must
be provided in order to make the steppers work. In the following we will
give some examples on reasonable state_typealgebraoperations combinations.
For the most common state types, like vector<double>
or array<double,N>
the default values range_algebra and default_operations are perfectly fine
and odeint can be used as is without worrying about algebra/operations at
all.
Important  

state_type, algebra and operations are not independent, a valid combination must be provided to make odeint work properly 
Moreover, as odeint handles the memory required for intermediate temporary objects itself, it also needs knowledge about how to create state_type objects and maybe how to allocate memory (resizing). All in all, the following things have to be taken care of when odeint is used with nonstandard state types:
Again, odeint already provides basic interfaces for most of the usual state
types. So if you use a std::vector
,
or a boost::array
as state type no additional work
is required, they just work out of the box.
We distinguish between two basic state types: fixed sized and dynamically
sized. For fixed size state types the default constructor state_type()
already allocates the required memory, prominent example is boost::array<T,N>
. Dynamically sized types have to be
resized to make sure enough memory is allocated, the standard constructor
does not take care of the resizing. Examples for this are the STL containers
like vector<double>
.
The most easy way of getting your own state type to work with odeint is to use a fixed size state, base calculations on the range_algebra and provide the following functionality:
Name 
Expression 
Type 
Semantics 

Construct State 


Creates an instance of 
Begin of the sequence 
boost::begin(x) 
Iterator 
Returns an iterator pointing to the begin of the sequence 
End of the sequence 
boost::end(x) 
Iterator 
Returns an iterator pointing to the end of the sequence 
Warning  

If your state type does not allocate memory by default construction, you must define it as resizeable and provide resize functionality (see below). Otherwise segmentation faults will occur. 
So fixed sized arrays supported by Boost.Range
immediately work with odeint. For dynamically sized arrays one has to additionally
supply the resize functionality. First, the state has to be tagged as resizeable
by specializing the struct is_resizeable
which consists of one typedef and one bool value:
Name 
Expression 
Type 
Semantics 

Resizability 


Determines resizeability of the state type, returns 
Resizability 


Same as above, but with 
Defining type
to be true_type
and value
as true
tells odeint that
your state is resizeable. By default, odeint now expects the support of
boost::size(x)
and
a x.resize( boost::size(y) )
member function for resizing:
Name 
Expression 
Type 
Semantics 

Get size 


Returns the current size of x. 
Resize 


Resizes x to have the same size as y. 
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 outofbox. 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< int 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
zerosized 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:
Name 
Expression 
Type 
Semantics 

Check size 


Returns true if the size of x equals the size of y. 
Resize 


Resizes x to have the same size as y. 
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::list
s
and so they can be used as state types. A complete example can be found
in list_lattice.cpp.
To provide maximum flexibility odeint is implemented in a highly modularized
way. This means it is possible to change the underlying mathematical operations
without touching the integration algorithms. The fundamental mathematical
operations are those of a vector space, that is addition of state_types
and multiplication of state_type
s with a scalar (time_type
). In odeint this is realized
in two concepts: Algebra and Operations. The standard way how this works
is by the range algebra which provides functions that apply a specific
operation to each of the individual elements of a container based on the
Boost.Range
library. If your state type is not supported by Boost.Range
there are several possibilities to tell odeint how to do algebraic operations:
boost::begin
and boost::end
for your state type so it works with Boost.Range.
+
and scalarvector multiplication operator *
and use the nonstandard vector_space_algebra
.
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_vector
s
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 
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_vector
s 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 precompiled
C library.
As seen above, the standard way of performing algebraic operations on
containerlike state types in odeint is to iterate through the elements
of the container and perform the operations elementwise 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:
Name 
Expression 
Type 
Semantics 

Addition 


Calculates the vector sum 'x+y'. 
Assign addition 


Performs x+y in place. 
Scalar multiplication 


Performs multiplication of vector x with scalar a. 
Assign scalar multiplication 


Performs inplace multiplication of vector x with scalar a. 
Defining these operators makes your state type work with any basic RungeKutta stepper. However, if you want to use stepsize 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:
Name 
Expression 
Type 
Semantics 

Division 


Calculates the elementwise division 'x/y' 
Absolute value 


Element wise absolute value 
Reduce 


Performs the

As an example for the employment of the vector_space_algebra
we will adopt ublas::vector
from Boost.uBLAS
to work as a state type in odeint. This is particularly easy because
ublas::vector
supports vectorvector addition
and scalarvector multiplication described above as well as boost::size
. It also has a resize member function
so all that has to be done in this case is to declare resizability:
typedef boost::numeric::ublas::vector< double > state_type; namespace boost { namespace numeric { namespace odeint { template<> struct is_resizeable< state_type > { typedef boost::true_type type; const static bool value = type::value; }; } } }
Now ublas::vector can be used as state type for simple RungeKutta steppers
in odeint by specifying the vector_space_algebra
as algebra in the template parameter list of the stepper. The following
code shows the corresponding definitions:
int main() { state_type x(3); x[0] = 10.0; x[1] = 5.0 ; x[2] = 0.0; typedef runge_kutta4< state_type , double , state_type , double , vector_space_algebra > stepper; integrate_const( stepper() , lorenz , x , 0.0 , 10.0 , 0.1 ); }
Note again, that we haven't supported the requirements for controlled steppers, but only for simple RungeKutta methods. You can find the full example in lorenz_ublas.cpp.
Here we show how to implement the required operators on a state type.
As example we define a new class point3D
representing a threedimensional 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 RungeKutta 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 elementwise 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/p1.z ); } point3D abs( const point3D &p ) { return point3D( std::abs(p.x) , std::abs(p.y) , std::abs(p.z) ); }
Finally, we have to add a specialization for reduce
implementing a reduction over the state type:
namespace boost { namespace numeric { namespace odeint { // specialization of vector_space_reduce, only required for steppers with error control template<> struct vector_space_reduce< point3D > { template< class Value , class Op > Value operator() ( const point3D &p , Op op , Value init ) { init = op( init , p.x ); //std::cout << init << " "; init = op( init , p.y ); //std::cout << init << " "; init = op( init , p.z ); //std::cout << init << std::endl; return init; } }; } } }
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 using the vector_space_algebra
again:
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>( 1E10 , 1E10 ) , 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
gsl_vector, gsl_matrix, ublas::matrix, blitz::matrix, thrust