Of course, odeint can handle complex state types, hence ODEs which are defined on complex vector spaces. An example is the Stuart-Landau oscillator

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

where Ψ and i is a complex variable. Such systems can easily be represent by an state type like array< complex< double > , 1 >. The definition of this ODE in C++ code is very simple

typedef boost::array< complex< double > , 1 > 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[0] = ( 1.0 + m_eta * I ) * x[0] - ( 1.0 + m_alpha * I ) * norm( x[0] ) * x[0];
    }
};

Of courst, 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.

Integration is also very easy:

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

const double dt = 0.1;

typedef runge_kutta4< state_type > 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

odeint can also be used to solved 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.

The FPU is solved again by a symplectic solver, but in our 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 peformance 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 std::tr1::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 is more or less trivial

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 oscialltor 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 strenght 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 type unit and dimension analysis. The main problem here is, that the state_type is now heterogeneous, meaning that every entry has a different type. To overcome this problem compile-type 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. They have to be explicitly defined. Next, we define the ordinary differential equation which is completety equivalent to the example in the first section of this tutorial

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 );
    }
};

Having done this, the most difficulties have been done. Next, we instantiate an appropriate stepper. We must explicitly parametrize the stepper with the state_type, deriv_type, time_type. Furthermore, the basic calculations are now performed by the fusion_algebra which must also be given.

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.

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 phas oscillators, they are extremly 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.

Wave equation

KdV

Ginzburg-Landau

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 carful 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 used

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 type'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 exemplarily show this for a Hamiltonian system of nonlinear, disordered oscillators with nonlinear nearest neighbour 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) );

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