Is boost odeint supporting dimensional analysis via Boost.Units?

Question

I was trying to implement a simple ode to test whether boost.odeint is supporting the usage of boost.units. However my example is failing at compilation. Is it my code, or doesn't boost.odeint support boost.Units for dimensional analysis?

#include<boost/units/systems/si.hpp>
#include<boost/numeric/odeint.hpp>

typedef boost::units::quantity<boost::units::si::length,double> state_type;
typedef boost::units::quantity<velocity,double> dev_type;
typedef boost::units::quantity<boost::units::si::time,double> time_type;


using namespace boost::units::si;

void exponential_decay(const state_type &x, dev_type &dxdt, time_type t){
    dxdt = -0.0001*x/second;
}

int main(){
    state_type startValue = 10*meter;

    using namespace boost::numeric::odeint;
    auto steps = integrate(exponential_decay, startValue, 0.0*second,
            10.0*second, 0.1*second);

    return 0;
}

Answer 1

So, after initially failing to find a working set of state_type, algebra and operations¹ @Arash supplied the missing link.

However, for completeness, you do not need to bring the heavy machinery (fusion_algebra.hpp). These are for the more general cases, e.g. where the state type is not a single value.

In this simple case, all that was really required was to specify the algebra, instead of going with the default. This involves declaring a stepper_type, like:

using stepper_type = runge_kutta4<length_type, double, velocity_type,
                            time_type, vector_space_algebra>;

Next, we pick an itegration algorithm overload that allows us to supply that stepper:

Live On Coliru

#include <boost/numeric/odeint.hpp>
#include <boost/units/systems/si.hpp>

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

namespace si = boost::units::si;

void exponential_decay(const length_type &x, velocity_type &dxdt, time_type /*t*/) { dxdt = -0.0001 * x / si::second; }

int main() {
    using stepper_type = boost::numeric::odeint::runge_kutta4<
        length_type, double, velocity_type, time_type, boost::numeric::odeint::vector_space_algebra>;

    length_type startValue = 10 * si::meter;
    auto steps = integrate_const(
            stepper_type{}, exponential_decay,
            startValue, 0.0 * si::second, 10.0 * si::second,
            0.1 * si::second);

    std::cout << "Steps: " << steps << "\n";
}

Prints

Steps: 100

---

¹ from just looking at http://www.boost.org/doc/libs/1_66_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html and http://www.boost.org/doc/libs/1_66_0/doc/html/boost_units/Quantities.html

Answer 2

Use boost::fusion to fix the problem.

#include <iostream>
#include <vector>

#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/algebra/fusion_algebra.hpp>
#include <boost/numeric/odeint/algebra/fusion_algebra_dispatcher.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  > state_type;
typedef fusion::vector< velocity_type > deriv_type;

void exponential_decay( const state_type &x , deriv_type &dxdt , time_type t )
{
    fusion::at_c< 0 >( dxdt ) = (-0.0001/si::second)* fusion::at_c< 0 >( x );
}

void observer(const state_type &x,const time_type &t )
{
}

int main( int argc , char**argv )
{
    typedef runge_kutta_dopri5< state_type , double , deriv_type , time_type > stepper_type;

    state_type startValue( 1.0 * si::meter );

    auto steps=integrate_adaptive(
        make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type()),
        exponential_decay,
        startValue , 0.0 * si::second , 100.0 * si::second , 0.1 * si::second , observer);

    std::cout<<"steps: "<<steps<<std::endl;

    return 0;
}