...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 |

Imaging, you want to numerically integrate a harmonic oscillator with friction.
The equations of motion are given by *x'' = -x + γ x'*.
Odeint only deals with first order ODEs that have no higher derivatives than
x' involved. However, any higher order ODE can be transformed to a system
of first order ODEs by introducing the new variables *q=x*
and *p=x'* such that *w=(q,p)*. To
apply numerical integration one first has to design the right hand side of
the equation *w' = f(w) = (p,-q+γ p)*:

/* The type of container used to hold the state vector */ typedef std::vector< double > state_type; const double gam = 0.15; /* The rhs of x' = f(x) */ void harmonic_oscillator( const state_type &x , state_type &dxdt , const double /* t */ ) { dxdt[0] = x[1]; dxdt[1] = -x[0] - gam*x[1]; }

Here we chose `vector<double>`

as the state type, but others are also possible, for example `boost::array<double,2>`

. odeint is designed in such a way that
you can easily use your own state types. Next, the ODE is defined which is
in this case a simple function calculating *f(x)*. The
parameter signature of this function is crucial: the integration methods
will always call them in the form ```
f(x,
dxdt,
t)
```

(there are exceptions for some special routines). So, even if there is no
explicit time dependence, one has to define `t`

as a function parameter.

Now, we have to define the initial state from which the integration should start:

state_type x(2); x[0] = 1.0; // start at x=1.0, p=0.0 x[1] = 0.0;

For the integration itself we'll use the `integrate`

function, which is a convenient way to get quick results. It is based on
the error-controlled `runge_kutta54_cash_karp`

stepper (5th order) and uses adaptive step-size.

size_t steps = integrate( harmonic_oscillator , x , 0.0 , 10.0 , 0.1 );

The integrate function expects as parameters the rhs of the ode as defined
above, the initial state `x`

,
the start-and end-time of the integration as well as the initial time step=size.
Note, that `integrate`

uses
an adaptive step-size during the integration steps so the time points will
not be equally spaced. The integration returns the number of steps that were
applied and updates x which is set to the approximate solution of the ODE
at the end of integration.

It is also possible to represent the ode system as a class. The rhs must then be implemented as a functor - a class with an overloaded function call operator:

/* The rhs of x' = f(x) defined as a class */ class harm_osc { double m_gam; public: harm_osc( double gam ) : m_gam(gam) { } void operator() ( const state_type &x , state_type &dxdt , const double /* t */ ) { dxdt[0] = x[1]; dxdt[1] = -x[0] - m_gam*x[1]; } };

which can be used via

harm_osc ho(0.15); steps = integrate( ho , x , 0.0 , 10.0 , 0.1 );

In order to observe the solution during the integration steps all you have to do is to provide a reasonable observer. An example is

struct push_back_state_and_time { std::vector< state_type >& m_states; std::vector< double >& m_times; push_back_state_and_time( std::vector< state_type > &states , std::vector< double > × ) : m_states( states ) , m_times( times ) { } void operator()( const state_type &x , double t ) { m_states.push_back( x ); m_times.push_back( t ); } };

which stores the intermediate steps in a container. Note, the argument structure of the ()-operator: odeint calls the observer exactly in this way, providing the current state and time. Now, you only have to pass this container to the integration function:

vector<state_type> x_vec; vector<double> times; steps = integrate( harmonic_oscillator , x , 0.0 , 10.0 , 0.1 , push_back_state_and_time( x_vec , times ) ); /* output */ for( size_t i=0; i<=steps; i++ ) { cout << times[i] << '\t' << x_vec[i][0] << '\t' << x_vec[i][1] << '\n'; }

That is all. You can use functional libraries like Boost.Lambda or Boost.Phoenix to ease the creation of observer functions.

The full cpp file for this example can be found here: harmonic_oscillator.cpp