Solving ordinary differential equations in C++

• Download source code - 40.57 KB

Attention: A new version of odeint exists, which is decribed here.

Introduction

This article introduces the C++ framework odeint for solving ordinary differential equations (ODEs), which is based on template meta-programming. I think this framework has some nice advantages over existing code on ODEs, and it uses templates in a very elegant way. Furthermore, odeint has a STL-like syntax, and is independent of a specific container. The whole project lives in the boost sandbox; the code presented here is more or less a snapshot of the current development.

Background

I will only very briefly describe ordinary differential equations. If there is some interest in a more detailed explanation of ODEs, I can extend this part in future versions of the article.

Ordinary Differential Equations

An ordinary differential equation describes the evolution of some quantity x in terms of its derivative. It often takes the form:

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

The function f defines the ODE, and x and f can be vectors. Associated with every ODE is an initial value problem (IVP) that is the ODE, and 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, and it can be shown that every IVP has a unique solution. Of course, ordinary differential equations are not restricted to temporal problems, hence the variable t can be replaced by another quantity, like a spatial coordinate. ODEs and their relative PDEs (partial differential equation) are very important in nearly all scientific disciplines. All major equations in physics fall in this class, like Newton's law for classical physics, the Maxwell's equations for electromagnetism, the Schrödinger equation and its relativistic generalizations for the quantum world, or Einstein's equation for general relativity. Very often, the spatial axes of a PDE are discretized, and an ODE on a lattice results.

In the following figure, an example of an ODE from chaos theory is shown: the famous Lorenz attractor. It was the first example of a deterministic chaotic system, and triggered a huge amount of scientific work. The Lorenz attractor is described by a set of coupled ordinary differential equations:

d x1 / dt = sigma * ( x2 - x1 )
d x2 / dt = R * x1 - x2 - x1 * x3
d x3 / dt = x1 * x2 - b * x3

with:

sigma = 10
R = 28
b = 8 / 3

This equation has no analytical solution, such that it can only be solved numerically.

Numerical Solutions of ODEs

To solve ODEs numerically, various methods exist; all of them discretize the time. The easiest method is surely the explicit Euler scheme, which writes the derivative as the difference quotient:

d x(t) / d t = x(t+dt) - x(t) / dt

Then, some basic algebraic manipulations yield:

x(t+dt) = x(t) + dt*f (x(t) , t )

such that x(t+dt) can be obtained from the knowledge of x at the time t. Note again, that x does not have to be a scalar, but can be a vector. The rest is very simple. One chooses an initial condition x(t=0)=x0, and iteratively applies the Euler step to x(t) to obtain the evolution of x(t). Of course, more advanced solvers exist, and the most commonly used solver is probably the Runge-Kutta method of fourth order.

In the following source code, the function f(x(t)) defining the ODE is synonymously called system or dynamical system.

The Main Ideas

Now, we come to the design aspects of the library. One clearly needs two main building blocks:

1. A stepper function, which calculates the state x(t+d t) from the knowledge of x(t).
2. An integrate function, which performs the iteration, hence applies the stepper function many times.

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

Stepper Methods

We try to avoid abstract classes, since this costs some extra performance. Instead, we use generic programming and concepts which describe how a type or class has to behave, and whose methods have to be available. For the stepper function, we need the following ingredients:

template<
    class Container ,
    class Time = double ,
    class Traits = container_traits< Container >
    >
class ode_step
{
    // provide basic typedefs
    //
 public:

    typedef unsigned short order_type;
    typedef Time time_type;
    typedef Traits traits_type;
    typedef typename traits_type::container_type container_type;
    typedef typename traits_type::value_type value_type;

    // public interface
 public:

    // return the order of the stepper
    order_type order_step() const;

    // standard constructor
    ode_step( void );

    // adjust the size of the internal arrays
    void adjust_size( const container_type &x );

    // performs one step
    template< class DynamicalSystem >
    void do_step( DynamicalSystem &system ,
                  container_type &x ,
                  const container_type &dxdt ,
                  time_type t ,
                  time_type dt );

    // performs one step
    template< class DynamicalSystem >
    void do_step( DynamicalSystem &system ,
                  container_type &x ,
                  time_type t ,
                  time_type dt );
};

The container type defines how the state is stored. It has to fulfill some basic requirements like the value type and it should be default constructible. The container_traits type abstracts the resize functionality of the sequence concept. Furthermore, it knows how to obtain the appropriate begin() and end() iterators. Its usage will be shown below. The adust_size method is responsible for adjusting the size of any internal containers, which store the results of intermediate calculations.

The do_step functions calculates x(t+dt) from the knowledge of x(t), and they transform x(t) in-place, meaning that x(t) is replaced by x(t+dt) within do_step. The argument system defines the ODE. This can be a simple function with three arguments: system( state , dxdt , t ), or a class with the appropriate brackets operator()(state,dxdt,t). Such a class is also known as a functor. Two versions of do_step exist: one which calculates x(t+dt) with the explicit knowledge of f(x(t)), and one which calculates this value internally. In principle, the second version only calls the first version. Both functions exist, because in many situations, the user needs the knowledge of f(x(t)) and will call the system method with x(t). To avoid a double call of system, the do_step version with the explicit declaration of f(x(t)) is introduced.

Note, that the above class definition is not found in the code; it is only used here to show the stepper functionality. All specific classes possessing these typedefs and methods are said to be stepper classes; we also say that all classes with these definitions fulfill the stepper concept.

Now, we look at a specific method, the Euler solver:

template<
    class Container ,
    class Time = double ,
    class Traits = container_traits< Container >
    >
class stepper_euler
{
    // provide basic typedefs
 public:

    typedef unsigned short order_type;
    typedef Time time_type;
    typedef Traits traits_type;
    typedef typename traits_type::container_type container_type;
    typedef typename traits_type::value_type value_type;

    // private members
 private:

    container_type m_dxdt;

    // public interface
 public:

    // return the order of the stepper
    order_type order_step() const { return 1; }

    // standard constructor, m_dxdt is not adjusted
    stepper_euler( void )
        {
        }

    // contructor, which adjusts m_dxdt
    stepper_euler( const container_type &x )
    {
        adjust_size( x );
    }

    // adjust the size of m_dxdt
    void adjust_size( const container_type &x )
    {
        traits_type::adjust_size( x , m_dxdt );
    }

    // performs one step with the knowledge of dxdt(t)
    template< class DynamicalSystem >
    void do_step( DynamicalSystem &system ,
                  container_type &x ,
                  const container_type &dxdt ,
                  time_type t ,
                  time_type dt )
    {
        detail::it_algebra::increment( traits_type::begin(x) ,
                                       traits_type::end(x) ,
                                       traits_type::begin(dxdt) ,
                                       dt );
    }

    // performs one step
    template< class DynamicalSystem >
    void do_step( DynamicalSystem &system ,
                  container_type &x ,
                  time_type t ,
                  time_type dt )
    {
        system( x , m_dxdt , t );
        do_step( system , x , m_dxdt , t , dt );
    }
};

This definition is straightforward, only some notes:

Adjusting the size is done in the container_traits-type. The standard traits type is defined as:

template< class Container >
struct container_traits
{
    typedef Container container_type;
    typedef typename container_type::value_type value_type;
    typedef typename container_type::iterator iterator;
    typedef typename container_type::const_iterator const_iterator;

    static void resize( const container_type &x , container_type &dxdt ) 
	{ dxdt.resize( x.size() ); }
    static bool same_size( const container_type &x1 , const container_type &x2 ) 
	{ return (x1.size() == x2.size()); }
    static void adjust_size( const container_type &x1 , container_type &x2 ) 
	{ if( !same_size( x1 , x2 ) ) resize( x1 , x2 ); }

    static iterator begin( container_type &x ) { return x.begin(); }
    static const_iterator begin( const container_type &x ) { return x.begin(); }
    static iterator end( container_type &x ) { return x.end(); }
    static const_iterator end( const container_type &x ) { return x.end(); }
};

which will work with all (STL) containers, fulfilling the sequence and the container concept. But for other containers, these requirements are not sufficient. For example std::tr1::array does not possess a resize method. The container traits are also responsible to for the begin() and end() iterators. For SGI containers this is trivial, but for matrix types like Blitz++ or MTL matrices the begin() and end() are not defined, but can be obtained in different ways.

All iterator operations have been put in their own functions. This has the advantage that the code is small and clear, since no iterators need to be defined in the stepper functions. Furthermore, it is be possible to specialize and optimize these iterator functions for specific iterators or containers.

The adjust size functionality is necessary if the system size changes during the evolution; think, for example, of an ODE defined on a network with a non-constant number of nodes.

This is not the end of the story of the stepper functions. odeint also provides an advanced stepper concept which calculates the numerical errors made during each step. With this error calculation, it is possible to adapt the step-size dt during each step in the iteration. The concepts behind the adaptive step size integration will not be described in this article.

The following code demonstrates the usage of the stepper function:

#include <iostream>
#include <boost/numeric/odeint.hpp>

#define tab "\t"

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

typedef vector< double > state_type;

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

void lorenz( 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] = x[0]*x[1] - b * x[2];
}

int main( int argc , char **argv )
{
    const double dt = 0.01;

    state_type x(3);
    x[0] = 1.0 ;
    x[1] = 0.0 ;
    x[2] = 0.0;
    stepper_euler< state_type > stepper;
    stepper.adjust_size( x );

    double t = 0.0;
    for( size_t oi=0 ; oi<10000 ; ++oi,t+=dt )
    {
        stepper.do_step( lorenz , x , t , dt );
        cout << x[0] << tab << x[1] << tab << x[2] << endl;
    }
}

In some cases, it is desirable to define the ODE as a class. In our example, lorenz could also be defined as a class, having the appropriate bracket operator:

class lorenz_class
{
 public:
    void operator()( 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] = x[0]*x[1] - b * x[2];
    }
}

which is then called via:

lorenz_class lorenzo;
stepper.do_step( lorenzo , x , t , dt );

With classes more complicated dynamical systems (ODEs) can be created.

Integration Function

The next main idea is to provide one or more functions to automatically iterate the state of the ODE, with the possibility to observe the state in each step. This is done with the integrate_const function:

template<
    class Stepper ,
    class DynamicalSystem ,
    class Observer
>
size_t integrate_const(
    Stepper &stepper ,
    DynamicalSystem &system ,
    typename Stepper::container_type &state ,
    typename Stepper::time_type start_time ,
    typename Stepper::time_type end_time ,
    typename Stepper::time_type dt ,
    Observer &observer
    )
{
    stepper.adjust_size( state );
    size_t iteration = 0;
    while( start_time < end_time )
    {
        observer( start_time , state , system );
        stepper.do_step( system , state , start_time , dt );
        start_time += dt;
        ++iteration;
    }
    observer( start_time , state , system );

    return iteration;
}

which solves the ODE with constant steps dt in the time interval start_time, end_time. The initial state has also to be provided. Furthermore, the concept of an observer is introduced, although this is not exactly an Observer pattern. Again, the observer can have any type, it only has to possesses the bracket operator, operator()( time , state , system ). For example, lambda expressions from boost can be used:

integrate_const( stepper , lorenz , x , 0.0 , 10.0 , dt ,
    cout << _1 << tab << _2[0] << "\n" );

or:

vector< double > lx;
back_insert_iterator< vector<double> > iter( lx );
integrate_const( stepper , lorenz , x , 0.0 , 10.0 , dt , var(*iter++) = _2[0] );

Installation

The whole library is header only, meaning that no special installation process has to be carried out. On unixoid systems, do the following:

1. Unpack the sources : unzip odeint.zip
2. Change to the odeint directory : cd odeint
3. Set the boost environment variable : export $BOOST_ROOT=path_to_boost
4. Compile the lorenz example : make

Now, an executable lorenz should be present. The boost sources are needed for the lambda expressions in the example. You can download them from boost.org and the environment variable should point to the root directory of boost, in which you will find the subdirectories bin.v2, boost, doc, lib, more, ...

For MSVC odeint will also work, but at the moment I don't have a step to step guide on how to compile the example.

Summary and Outlook

I have shown you some design aspects, and the usage of odeint - a C++ framework for solving ordinary differential equations. It is very easy to use and very flexible. The performance is also very good. We have made some test run, where odeint is compared with the gsl. As a test system, the Lorenz system was used. It turned out that odeint is about two times faster than the gsl routines.

The source code provided with this article is a snapshot of the current development. Some methods for adaptive step size control exist and some methods to solve Hamiltonian system. In the following table, an overview over existing methods is given. S means that the stepper fulfills the simple stepper concept described above, whereas E means that the stepper fulfills the Error stepper concept, which is not introduced here but needed for adaptive step size control.

Method	Class name	Order	Error (order)	Stepper concept
Euler	stepper_euler	1	No	S
Runge Kutta 4	stepper_rk4	4	No	S
Runge Kutta Cash Karp	stepper_rk5_ck	5	Yes 4	E
Runge Kutta Fehlber	stepper_rk78_fehlberg	7	Yes 8	S,E
Midpoint	stepper_midpoint	Var.	No	S
Symplectic Euler	hamiltonian_stepper_euler	1	No	S
Symplectic Runge-Kutta	hamiltonian_stepper_rk	6	No	S

The whole library is not complete, and will be extended in the near future. We have a small Todo list, which should be completed before the final release:

• Adaptors for Blitz++, MTL, boost::ublas and boost::multi_array
• Documentation
• Solvers for stiff systems
• Abstract the iterator algebra

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

History

• 9.11.2009 - Initial document
• 16.3.2010 - Update, introducing the container traits, replacing the examples
• 23.3.2010 - Updated source code and article