13,355,726 members (63,230 online)
alternative version

#### Stats

192.4K views
96 bookmarked
Posted 10 Nov 2009

# Solving ordinary differential equations in C++

, 19 Oct 2011
 Rate this:
This article explains a framework for solving ordinary differential equations, which is based on template metaprogramming.

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 `typedef`s 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 )
{
}

stepper_euler( const container_type &x )
{
}

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

// 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;

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
)
{
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

## Share

 Germany
No Biography provided

## You may also be interested in...

 First PrevNext
 solve state-space model? Mohammed Alloboany30-Dec-15 1:47 Mohammed Alloboany 30-Dec-15 1:47
 My vote of 5 Edward Keningham11-Aug-12 4:22 Edward Keningham 11-Aug-12 4:22
 My vote of 1 RedDK22-Feb-12 7:30 RedDK 22-Feb-12 7:30
 Excellent stuff Espen Harlinn11-Feb-12 6:07 Espen Harlinn 11-Feb-12 6:07
 speed question Antediluvian999923-Jan-12 6:27 Antediluvian9999 23-Jan-12 6:27
 Re: speed question Antediluvian999923-Jan-12 6:28 Antediluvian9999 23-Jan-12 6:28
 Re: speed question Antediluvian999923-Jan-12 6:47 Antediluvian9999 23-Jan-12 6:47
 Re: speed question Antediluvian999923-Jan-12 7:02 Antediluvian9999 23-Jan-12 7:02
 Re: speed question Antediluvian999915-Feb-12 15:58 Antediluvian9999 15-Feb-12 15:58
 Re: speed question Antediluvian999917-Feb-12 0:46 Antediluvian9999 17-Feb-12 0:46
 My vote of 1 RedDK20-Oct-11 6:55 RedDK 20-Oct-11 6:55
 Re: My vote of 1 AndyUk0624-Oct-11 5:33 AndyUk06 24-Oct-11 5:33
 Re: My vote of 1 RedDK24-Oct-11 7:04 RedDK 24-Oct-11 7:04
 Adaptive Example safariman11-May-10 8:55 safariman 11-May-10 8:55
 Re: Adaptive Example safariman13-May-10 7:24 safariman 13-May-10 7:24
 Thank you. The code snippet you provided works. However, when I change `typedef stepper_rk5_ck< state_type > stepper_type;` to `typedef stepper_rk78_fehlberg< state_type > stepper_type;` I get the long-winded error pasted below that basically says that rk78_fehlberg does not have function `do_step()`. Is this a bug with my code or is controlled stepping not implemented in rk78_fehlberg yet? Thanks. ```/Users/philipadmin/Packages/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp: In member function 'boost::numeric::odeint::controlled_step_result boost::numeric::odeint::controlled_stepper_standard::try_step(DynamicalSystem&, typename ErrorStepper::container_type&, const typename ErrorStepper::container_type&, typename ErrorStepper::time_type&, typename ErrorStepper::time_type&) [with DynamicalSystem = VehicleDynamics, ErrorStepper = boost::numeric::odeint::stepper_rk78_fehlberg >, double, boost::numeric::odeint::container_traits > > >]': /Users/philipadmin/Packages/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp:199: instantiated from 'boost::numeric::odeint::controlled_step_result boost::numeric::odeint::controlled_stepper_standard::try_step(DynamicalSystem&, typename ErrorStepper::container_type&, typename ErrorStepper::time_type&, typename ErrorStepper::time_type&) [with DynamicalSystem = VehicleDynamics, ErrorStepper = boost::numeric::odeint::stepper_rk78_fehlberg >, double, boost::numeric::odeint::container_traits > > >]' /Users/philipadmin/Packages/odeint/boost/numeric/odeint/integrator_adaptive_stepsize.hpp:54: instantiated from 'size_t boost::numeric::odeint::integrate_adaptive(ControlledStepper&, DynamicalSystem&, typename ControlledStepper::container_type&, typename ControlledStepper::time_type, typename ControlledStepper::time_type, typename ControlledStepper::time_type, Observer&) [with ControlledStepper = boost::numeric::odeint::controlled_stepper_standard >, double, boost::numeric::odeint::container_traits > > > >, DynamicalSystem = VehicleDynamics, Observer = void ()(double, std::vector >&, VehicleDynamics&)]' /Users/philipadmin/Packages/odeint/boost/numeric/odeint/integrator_adaptive_stepsize.hpp:87: instantiated from 'size_t boost::numeric::odeint::integrate_adaptive(ControlledStepper&, DynamicalSystem&, typename ControlledStepper::container_type&, typename ControlledStepper::time_type, typename ControlledStepper::time_type, typename ControlledStepper::time_type) [with ControlledStepper = boost::numeric::odeint::controlled_stepper_standard >, double, boost::numeric::odeint::container_traits > > > >, DynamicalSystem = VehicleDynamics]' ../src/TrueVehicle.cpp:34: instantiated from here /Users/philipadmin/Packages/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp:158: error: no matching function for call to 'boost::numeric::odeint::stepper_rk78_fehlberg >, double, boost::numeric::odeint::container_traits > > >::do_step(VehicleDynamics&, std::vector >&, const std::vector >&, double&, double&, std::vector >&)' make: *** [src/TrueVehicle.o] Error 1```