Thus far we have seen several examples defined for real values. Of course, odeint can handle complex state types, hence ODEs which are defined on complex vector spaces, as well. An example is the Stuart-Landau oscillator

d Ψ / dt = ( 1 + i η ) Ψ + ( 1 + i α ) | Ψ |² Ψ

where Ψ and i is a complex variable. Hence the state type of this system is complex< double >. The definition of this ODE in C++ code is very simple

typedef complex< double > state_type;

struct stuart_landau
{
    double m_eta;
    double m_alpha;

    stuart_landau( double eta = 1.0 , double alpha = 1.0 )
    : m_eta( eta ) , m_alpha( alpha ) { }

    void operator()( const state_type &x , state_type &dxdt , double t ) const
    {
        const complex< double > I( 0.0 , 1.0 );
        dxdt = ( 1.0 + m_eta * I ) * x - ( 1.0 + m_alpha * I ) * norm( x ) * x;
    }
};

Of course, one can also use classical functions to implement the state function. In this cast the Stuart-Landau oscillator looks like

double eta = 1.0;
double alpha = 1.0;

void stuart_landau( const state_type &x , state_type &dxdt , double t )
{
    const complex< double > I( 0.0 , 1.0 );
    dxdt[0] = ( 1.0 + m_eta * I ) * x[0] - ( 1.0 + m_alpha * I ) * norm( x[0] ) * x[0];
}

We strongly recommend to use the first ansatz. In this case you have explicit control over the parameters of the system and are not restricted to use global variables to parametrize the oscillator.

When chosing the stepper type one has to account for the "unusual" state type: it is a single complex<double> opposed to the vector types used in the previous examples. This means that no iterations over vector elements have to be performed inside the stepper algorithm. You can enforce this by supplying additional template arguments to the stepper including the vector_space_algebra. Details on the usage of algebras can be found in the section Adapt you own state types.

state_type x = complex< double >( 1.0 , 0.0 );

const double dt = 0.1;

typedef runge_kutta4< state_type , double , state_type , double ,
                      vector_space_algebra > stepper_type;

integrate_const( stepper_type() , stuart_landau( 2.0 , 1.0 ) , x , 0.0 , 10.0 , dt , streaming_observer( cout ) );

The full cpp file for the Stuart Landau example can be found here stuart_landau.cpp

	Note
	The fact that we have to configure a different algebra is solely due to the fact that we use a non-vector state type and not to the usage of complex values. So for, e.g. vector< complex<double> >, this would not be required.

odeint can also be used to solve ordinary differential equations defined on lattices. A prominent example is the Fermi-Pasta-Ulam system [8]. It is a Hamiltonian system of nonlinear coupled harmonic oscillators. The Hamiltonian is

H = Σ_​i p_​i²/2 + 1/2 ( q_​i+1 - q_​i )^2 + β / 4 ( q_​i+1 - q_​i )^4

Remarkably, the Fermi-Pasta-Ulam system was the first numerical experiment which has been implemented on a computer in 1953. It was studied at Los Alamos on one of the first computer (a MANIAC I) and it triggered a whole new tree of mathematical and physical science.

Like the Solar System, the FPU is solved again by a symplectic solver, but in this case we can speed up the computation because the q components trivially reduce to dq_​i / dt = p_​i. odeint is capable of doing this performance improvement. All you have to do is to call the symplectic solver with an state function for the p components. Here is how this function looks like

typedef vector< double > container_type;

struct fpu
{
    const double m_beta;

    fpu( const double beta = 1.0 ) : m_beta( beta ) { }

    // system function defining the ODE
    void operator()( const container_type &q , container_type &dpdt ) const
    {
        size_t n = q.size();
        double tmp = q[0] - 0.0;
        double tmp2 = tmp + m_beta * tmp * tmp * tmp;
        dpdt[0] = -tmp2;
        for( size_t i=0 ; i<n-1 ; ++i )
        {
            tmp = q[i+1] - q[i];
            tmp2 = tmp + m_beta * tmp * tmp * tmp;
            dpdt[i] += tmp2;
            dpdt[i+1] = -tmp2;
        }
        tmp = - q[n-1];
        tmp2 = tmp + m_beta * tmp * tmp * tmp;
        dpdt[n-1] += tmp2;
    }

    // calculates the energy of the system
    double energy( const container_type &q , const container_type &p ) const
    {
        // ...
    }

    // calculates the local energy of the system
    void local_energy( const container_type &q , const container_type &p , container_type &e ) const
    {
        // ...
    }
};

Of course, you can also use boost::array< double , N > for the state type.

Now, you have to define your initial values and perform the integration. All this can be easily done with the following piece of code:

const size_t n = 64;
container_type q( n , 0.0 ) , p( n , 0.0 );

for( size_t i=0 ; i<n ; ++i )
{
    p[i] = 0.0;
    q[i] = 32.0 * sin( double( i + 1 ) / double( n + 1 ) * M_PI );
}


const double dt = 0.1;

typedef symplectic_rkn_sb3a_mclachlan< container_type > stepper_type;
fpu fpu_instance( 8.0 );

integrate_const( stepper_type() , fpu_instance ,
        make_pair( boost::ref( q ) , boost::ref( p ) ) ,
        0.0 , 1000.0 , dt , streaming_observer( cout , fpu_instance , 10 ) );

The observer uses a reference to the system object to calculate the local energies:

struct streaming_observer
{
    std::ostream& m_out;
    const fpu &m_fpu;
    size_t m_write_every;
    size_t m_count;

    streaming_observer( std::ostream &out , const fpu &f , size_t write_every = 100 )
    : m_out( out ) , m_fpu( f ) , m_write_every( write_every ) , m_count( 0 ) { }

    template< class State >
    void operator()( const State &x , double t )
    {
        if( ( m_count % m_write_every ) == 0 )
        {
            container_type &q = x.first;
            container_type &p = x.second;
            container_type energy( q.size() );
            m_fpu.local_energy( q , p , energy );
            for( size_t i=0 ; i<q.size() ; ++i )
            {
                m_out << t << "\t" << i << "\t" << q[i] << "\t" << p[i] << "\t" << energy[i] << "\n";
            }
            m_out << "\n";
            clog << t << "\t" << accumulate( energy.begin() , energy.end() , 0.0 ) << "\n";
        }
        ++m_count;
    }
};

The full cpp file for this FPU example can be found here fpu.cpp

Another import high dimensional system of coupled ordinary differential equations is an ensemble of N all-to-all coupled phase oscillators [9]. It is defined as

dφ_​k / dt = ω_​k + ε / N Σ_​j sin( φ_​j - φ_​k )

The natural frequencies ω_​i of each oscillator follow some distribution and ε is the coupling strength. We choose here a Lorentzian distribution for ω_​i. Interestingly a phase transition can be observed if the coupling strength exceeds a critical value. Above this value synchronization sets in and some of the oscillators oscillate with the same frequency despite their different natural frequencies. The transition is also called Kuramoto transition. Its behavior can be analyzed by employing the mean field of the phase

Z = K e^{i Θ} = 1 / N Σ_​ke^{i φ_​k}

The definition of the system function is now a bit more complex since we also need to store the individual frequencies of each oscillator.

typedef vector< double > container_type;


pair< double , double > calc_mean_field( const container_type &x )
{
    size_t n = x.size();
    double cos_sum = 0.0 , sin_sum = 0.0;
    for( size_t i=0 ; i<n ; ++i )
    {
        cos_sum += cos( x[i] );
        sin_sum += sin( x[i] );
    }
    cos_sum /= double( n );
    sin_sum /= double( n );

    double K = sqrt( cos_sum * cos_sum + sin_sum * sin_sum );
    double Theta = atan2( sin_sum , cos_sum );

    return make_pair( K , Theta );
}


struct phase_ensemble
{
    container_type m_omega;
    double m_epsilon;

    phase_ensemble( const size_t n , double g = 1.0 , double epsilon = 1.0 )
    : m_omega( n , 0.0 ) , m_epsilon( epsilon )
    {
        create_frequencies( g );
    }

    void create_frequencies( double g )
    {
        boost::mt19937 rng;
        boost::cauchy_distribution<> cauchy( 0.0 , g );
        boost::variate_generator< boost::mt19937&, boost::cauchy_distribution<> > gen( rng , cauchy );
        generate( m_omega.begin() , m_omega.end() , gen );
    }

    void set_epsilon( double epsilon ) { m_epsilon = epsilon; }

    double get_epsilon( void ) const { return m_epsilon; }

    void operator()( const container_type &x , container_type &dxdt , double /* t */ ) const
    {
        pair< double , double > mean = calc_mean_field( x );
        for( size_t i=0 ; i<x.size() ; ++i )
            dxdt[i] = m_omega[i] + m_epsilon * mean.first * sin( mean.second - x[i] );
    }
};

Note, that we have used Z to simplify the equations of motion. Next, we create an observer which computes the value of Z and we record Z for different values of ε.

struct statistics_observer
{
    double m_K_mean;
    size_t m_count;

    statistics_observer( void )
    : m_K_mean( 0.0 ) , m_count( 0 ) { }

    template< class State >
    void operator()( const State &x , double t )
    {
        pair< double , double > mean = calc_mean_field( x );
        m_K_mean += mean.first;
        ++m_count;
    }

    double get_K_mean( void ) const { return ( m_count != 0 ) ? m_K_mean / double( m_count ) : 0.0 ; }

    void reset( void ) { m_K_mean = 0.0; m_count = 0; }
};

Now, we do several integrations for different values of ε and record Z. The result nicely confirms the analytical result of the phase transition, i.e. in our example the standard deviation of the Lorentzian is 1 such that the transition will be observed at ε = 2.

const size_t n = 16384;
const double dt = 0.1;

container_type x( n );

boost::mt19937 rng;
boost::uniform_real<> unif( 0.0 , 2.0 * M_PI );
boost::variate_generator< boost::mt19937&, boost::uniform_real<> > gen( rng , unif );

// gamma = 1, the phase transition occurs at epsilon = 2
phase_ensemble ensemble( n , 1.0 );
statistics_observer obs;

for( double epsilon = 0.0 ; epsilon < 5.0 ; epsilon += 0.1 )
{
    ensemble.set_epsilon( epsilon );
    obs.reset();

    // start with random initial conditions
    generate( x.begin() , x.end() , gen );

    // calculate some transients steps
    integrate_const( runge_kutta4< container_type >() , boost::ref( ensemble ) , x , 0.0 , 10.0 , dt );

    // integrate and compute the statistics
    integrate_const( runge_kutta4< container_type >() , boost::ref( ensemble ) , x , 0.0 , 100.0 , dt , boost::ref( obs ) );
    cout << epsilon << "\t" << obs.get_K_mean() << endl;
}

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

odeint also works well with Boost.Units - a library for compile time unit and dimension analysis. Ot works by decoding unit information into the types of values. For a one-dimensional unit you can just use the Boost.Unit types as state type, deriv type and time type and hand the vector_space_algebra to the stepper definition and everything works just fine:

typedef units::quantity< si::time , double > time_type;
typedef units::quantity< si::length , double > length_type;
typedef units::quantity< si::velocity , double > velocity_type;

typedef runge_kutta4< length_type , double , velocity_type , time_type ,
                      vector_space_algebra > stepper_type;

If you want to solve more-dimensional problems the individual entries typically have different units. That means that the state_type is now possibly heterogeneous, meaning that every entry might have a different type. To solve this problem, compile-time sequences from Boost.Fusion can be used.

To illustrate how odeint works with Boost.Units we use the harmonic oscillator as primary example. We start with defining all quantities

#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/algebra/fusion_algebra.hpp>

#include <boost/units/systems/si/length.hpp>
#include <boost/units/systems/si/time.hpp>
#include <boost/units/systems/si/velocity.hpp>
#include <boost/units/systems/si/acceleration.hpp>
#include <boost/units/systems/si/io.hpp>

#include <boost/fusion/container.hpp>

using namespace std;
using namespace boost::numeric::odeint;
namespace fusion = boost::fusion;
namespace units = boost::units;
namespace si = boost::units::si;

typedef units::quantity< si::time , double > time_type;
typedef units::quantity< si::length , double > length_type;
typedef units::quantity< si::velocity , double > velocity_type;
typedef units::quantity< si::acceleration , double > acceleration_type;
typedef units::quantity< si::frequency , double > frequency_type;

typedef fusion::vector< length_type , velocity_type > state_type;
typedef fusion::vector< velocity_type , acceleration_type > deriv_type;

Note, that the state_type and the deriv_type are now a compile-time fusion sequences. deriv_type represents x' and is now different from the state type as it has different unit definitions. Next, we define the ordinary differential equation which is completely equivalent to the example in Harmonic Oscillator:

struct oscillator
{
    frequency_type m_omega;

    oscillator( const frequency_type &omega = 1.0 * si::hertz ) : m_omega( omega ) { }

    void operator()( const state_type &x , deriv_type &dxdt , time_type t ) const
    {
        fusion::at_c< 0 >( dxdt ) = fusion::at_c< 1 >( x );
        fusion::at_c< 1 >( dxdt ) = - m_omega * m_omega * fusion::at_c< 0 >( x );
    }
};

Next, we instantiate an appropriate stepper. We must explicitly parametrize the stepper with the state_type, deriv_type, time_type. Furthermore, the iteration over vector elements is now done by the fusion_algebra which must also be given. For more on the state types / algebras see chapter Adapt you own state types.

typedef runge_kutta_dopri5< state_type , double , deriv_type , time_type , fusion_algebra > stepper_type;

state_type x( 1.0 * si::meter , 0.0 * si::meter_per_second );

integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() ) , oscillator( 2.0 * si::hertz ) , x , 0.0 * si::second , 100.0 * si::second , 0.1 * si::second , streaming_observer( cout ) );

It is quite easy but the compilation time might take very long. Furthermore, the observer is defined a bit different

struct streaming_observer
{
    std::ostream& m_out;

    streaming_observer( std::ostream &out ) : m_out( out ) { }

    struct write_element
    {
        std::ostream &m_out;
        write_element( std::ostream &out ) : m_out( out ) { };

        template< class T >
        void operator()( const T &t ) const
        {
            m_out << "\t" << t;
        }
    };

    template< class State , class Time >
    void operator()( const State &x , const Time &t ) const
    {
        m_out << t;
        fusion::for_each( x , write_element( m_out ) );
        m_out << "\n";
    }
};

	Caution
	Using Boost.Units works nicely but compilation can be very time and memory consuming. For example the unit test for Boost.Units take up to 4 GB of memory at compilation.

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

odeint works well with a variety of different state types. It is not restricted to pure vector-wise types, like vector< double >, array< double , N >, fusion::vector< double , double >, etc. but also works with types having a different topology then simple vectors. Here, we show how odeint can be used with matrices as states type, in the next section we will show how can be used to solve ODEs defined on complex networks.

By default, odeint can be used with ublas::matrix< T > as state type for matrices. A simple example is a two-dimensional lattice of coupled phase oscillators. We like phase oscillators, they are extremely easy and might serve for different demonstration purposes. Other matrix types like mtl::dense_matrix or blitz arrays and matrices can used as well but need some kind of activation in order to work with odeint. This activation is described in following sections,

The definition of the system is

typedef boost::numeric::ublas::matrix< double > state_type;

struct two_dimensional_phase_lattice
{
    two_dimensional_phase_lattice( double gamma = 0.5 )
    : m_gamma( gamma ) { }

    void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
    {
        size_t size1 = x.size1() , size2 = x.size2();

        for( size_t i=1 ; i<size1-1 ; ++i )
        {
            for( size_t j=1 ; j<size2-1 ; ++j )
            {
                dxdt( i , j ) =
                        coupling_func( x( i + 1 , j ) - x( i , j ) ) +
                        coupling_func( x( i - 1 , j ) - x( i , j ) ) +
                        coupling_func( x( i , j + 1 ) - x( i , j ) ) +
                        coupling_func( x( i , j - 1 ) - x( i , j ) );
            }
        }

        for( size_t i=0 ; i<x.size1() ; ++i ) dxdt( i , 0 ) = dxdt( i , x.size2() -1 ) = 0.0;
        for( size_t j=0 ; j<x.size2() ; ++j ) dxdt( 0 , j ) = dxdt( 0 , x.size1() -1 ) = 0.0;
    }

    double coupling_func( double x ) const
    {
        return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
    }

    double m_gamma;
};

This is in principle all. Please note, that the above code is far from being optimal. Better performance can be achieved if every interaction is only calculated once and iterators for columns and rows are used. Below are some visualizations of the evolution of this lattice equation.

The full cpp for this example can be found here two_dimensional_phase_lattice.cpp.

Besides the classical floating point number like float, double, complex< double > you can also use arbitrary precision types, like the types from gmp and mpfr. But you have to be careful about instantiating any numbers.

For gmp types you have to set the default precision before any number is instantiated. This can be easily done by calling mpf_set_default_prec( precision ) as the first function in your main program. Secondly, you can not use any global constant variables since they will not be set with the default precision you have already set.

Here is a simple example:

typedef mpf_class value_type;
typedef boost::array< value_type , 3 > state_type;

struct lorenz
{
    void operator()( const state_type &x , state_type &dxdt , value_type t ) const
    {
        const value_type sigma( 10.0 );
        const value_type R( 28.0 );
        const value_type b( value_type( 8.0 ) / value_type( 3.0 ) );

        dxdt[0] = sigma * ( x[1] - x[0] );
        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
        dxdt[2] = -b * x[2] + x[0] * x[1];
    }
};

which can be integrated:

const int precision = 1024;
mpf_set_default_prec( precision );

state_type x = {{ value_type( 10.0 ) , value_type( 10.0 ) , value_type( 10.0 ) }};

cout.precision( 1000 );
integrate_const( runge_kutta4< state_type , value_type >() ,
        lorenz() , x , value_type( 0.0 ) , value_type( 10.0 ) , value_type( value_type( 1.0 ) / value_type( 10.0 ) ) ,
        streaming_observer( cout ) );

	Caution
	The full support of arbitrary precision types depends on the functionality they provide. For example, the types from gmp are lacking of functions for calculating the power and arbitrary roots, hence they can not be used with the controlled steppers. In detail, for full support the min( x , y ), max( x , y ), pow( y , y ) must be callable.

The full example can be found at lorenz_gmpxx.cpp.

odeint supports changes of the state size during integration if a state_type is used which can be resized, like std::vector. The adjustment of the state's size has to be done from outside and the stepper has to be instantiated with always_resizer as the template argument for the resizer_type. In this configuration, the stepper checks for changes in the state size and adjust it's internal storage accordingly.

We exemplary show this for a Hamiltonian system of nonlinear, disordered oscillators with nonlinear nearest neighbor coupling.

The system function is implemented in terms of a class that also provides functions for calculating the energy. Note, that this class stores the random potential internally which is not resized, but rather a start index is kept which should be changed whenever the states' size change.

typedef vector< double > coord_type;
typedef pair< coord_type , coord_type > state_type;

struct compacton_lattice
{
    const int m_max_N;
    const double m_beta;
    int m_pot_start_index;
    vector< double > m_pot;

    compacton_lattice( int max_N , double beta , int pot_start_index )
        : m_max_N( max_N ) , m_beta( beta ) , m_pot_start_index( pot_start_index ) , m_pot( max_N )
    {
        srand( time( NULL ) );
        // fill random potential with iid values from [0,1]
        boost::mt19937 rng;
        boost::uniform_real<> unif( 0.0 , 1.0 );
        boost::variate_generator< boost::mt19937&, boost::uniform_real<> > gen( rng , unif );
        generate( m_pot.begin() , m_pot.end() , gen );
    }

    void operator()( const coord_type &q , coord_type &dpdt )
    {
        // calculate dpdt = -dH/dq of this hamiltonian system
        // dp_i/dt = - V_i * q_i^3 - beta*(q_i - q_{i-1})^3 + beta*(q_{i+1} - q_i)^3
        const int N = q.size();
        double diff = q[0] - q[N-1];
        for( int i=0 ; i<N ; ++i )
        {
            dpdt[i] = - m_pot[m_pot_start_index+i] * q[i]*q[i]*q[i] -
                    m_beta * diff*diff*diff;
            diff = q[(i+1) % N] - q[i];
            dpdt[i] += m_beta * diff*diff*diff;
        }
    }

    void energy_distribution( const coord_type &q , const coord_type &p , coord_type &energies )
    {
        // computes the energy per lattice site normalized by total energy
        const size_t N = q.size();
        double en = 0.0;
        for( size_t i=0 ; i<N ; i++ )
        {
            const double diff = q[(i+1) % N] - q[i];
            energies[i] = p[i]*p[i]/2.0
                + m_pot[m_pot_start_index+i]*q[i]*q[i]*q[i]*q[i]/4.0
                + m_beta/4.0 * diff*diff*diff*diff;
            en += energies[i];
        }
        en = 1.0/en;
        for( size_t i=0 ; i<N ; i++ )
        {
            energies[i] *= en;
        }
    }

    double energy( const coord_type &q , const coord_type &p )
    {
        // calculates the total energy of the excitation
        const size_t N = q.size();
        double en = 0.0;
        for( size_t i=0 ; i<N ; i++ )
        {
            const double diff = q[(i+1) % N] - q[i];
            en += p[i]*p[i]/2.0
                + m_pot[m_pot_start_index+i]*q[i]*q[i]*q[i]*q[i] / 4.0
                + m_beta/4.0 * diff*diff*diff*diff;
        }
        return en;
    }

    void change_pot_start( const int delta )
    {
        m_pot_start_index += delta;
    }
};

The total size we allow is 1024 and we start with an initial state size of 60.

//start with 60 sites
const int N_start = 60;
coord_type q( N_start , 0.0 );
q.reserve( max_N );
coord_type p( N_start , 0.0 );
p.reserve( max_N );
// start with uniform momentum distribution over 20 sites
fill( p.begin()+20 , p.end()-20 , 1.0/sqrt(20.0) );

coord_type distr( N_start , 0.0 );
distr.reserve( max_N );

// create the system
compacton_lattice lattice( max_N , beta , (max_N-N_start)/2 );

//create the stepper, note that we use an always_resizer because state size might change during steps
typedef symplectic_rkn_sb3a_mclachlan< coord_type , coord_type , double , coord_type , coord_type , double ,
        range_algebra , default_operations , always_resizer > hamiltonian_stepper;
hamiltonian_stepper stepper;
hamiltonian_stepper::state_type state = make_pair( q , p );

The lattice gets resized whenever the energy distribution comes close to the borders distr[10] > 1E-150, distr[distr.size()-10] > 1E-150. If we increase to the left, q and p have to be rotated because their resize function always appends at the end. Additionally, the start index of the potential changes in this case.

double t = 0.0;
const double dt = 0.1;
const int steps = 10000;
for( int step = 0 ; step < steps ; ++step )
{
    stepper.do_step( boost::ref(lattice) , state , t , dt );
    lattice.energy_distribution( state.first , state.second , distr );
    if( distr[10] > 1E-150 )
    {
        do_resize( state.first , state.second , distr , state.first.size()+20 );
        rotate( state.first.begin() , state.first.end()-20 , state.first.end() );
        rotate( state.second.begin() , state.second.end()-20 , state.second.end() );
        lattice.change_pot_start( -20 );
        cout << t << ": resized left to " << distr.size() << ", energy = " << lattice.energy( state.first , state.second ) << endl;
    }
    if( distr[distr.size()-10] > 1E-150 )
    {
        do_resize( state.first , state.second , distr , state.first.size()+20 );
        cout << t << ": resized right to " << distr.size() << ", energy = " << lattice.energy( state.first , state.second ) << endl;
    }
    t += dt;
}

The do_resize function simply calls vector.resize of q , p and distr.

void do_resize( coord_type &q , coord_type &p , coord_type &distr , const int N )
{
    q.resize( N );
    p.resize( N );
    distr.resize( N );
}

The full example can be found in resizing_lattice.cpp