# Boost C++ Libraries

...one of the most highly regarded and expertly designed C++ library projects in the world. Home Libraries People FAQ More

### Short Example

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 = x;
dxdt = -x - gam*x;
}
```

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 = 1.0; // start at x=1.0, p=0.0
x = 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 = x;
dxdt = -x - m_gam*x;
}
};
```

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 > &times )
: 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] << '\t' << x_vec[i] << '\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