Using Boost::odeint with Eigen::Matrix as state vector

Question

I'm trying to utilize the ODE integration capabilities of Boost using the Matrix class from Eigen 3 as my state vector, but I'm running into problems deep into Boost that I don't understand how to address.

A minimal example of what I'm trying to do:

#include <Eigen/Core>
#include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
#include <iostream>

using namespace Eigen;
using namespace boost::numeric::odeint;

template<size_t N>
using vector = Matrix<double, N, 1>;

typedef vector<3> state;

int main() {

    state X0;
    X0 << 1., 2., 3.;
    state xout = X0;

    runge_kutta_dopri5<state> stepper;

    // If I remove these lines, everything compiles fine
    stepper.do_step([](const state x, state dxdt, const double t) -> void { 
        dxdt = x;
    }, X0, 0.0, xout, 0.01);

    std::cout << xout << std::endl;
}

If I coment out the call to stepper.do_step everything compiles and runs just fine, but of course doesn't do anything interesting. If I don't, Boost vomits compile errors over my terminal, the first of which is

In file included from /usr/include/boost/mpl/aux_/begin_end_impl.hpp:20:0,
                 from /usr/include/boost/mpl/begin_end.hpp:18,
                 from /usr/include/boost/mpl/is_sequence.hpp:19,
                 from /usr/include/boost/fusion/support/detail/is_mpl_sequence.hpp:12,
                 from /usr/include/boost/fusion/support/tag_of.hpp:13,
                 from /usr/include/boost/fusion/support/is_sequence.hpp:11,
                 from /usr/include/boost/fusion/sequence/intrinsic_fwd.hpp:12,
                 from /usr/include/boost/fusion/sequence/intrinsic/front.hpp:10,
                 from /usr/include/boost/fusion/include/front.hpp:10,
                 from /usr/include/boost/numeric/odeint/util/is_resizeable.hpp:26,
                 from /usr/include/boost/numeric/odeint/util/state_wrapper.hpp:25,
                 from /usr/include/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp:27,
                 from /usr/include/boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp:24,
                 from /home/tlycken/exjobb/Code/alpha-orbit-follower/test/algebra/algebra-tests.cpp:2:
/usr/include/boost/mpl/eval_if.hpp: In instantiation of ‘struct boost::mpl::eval_if_c<true, boost::range_const_iterator<Eigen::Matrix<double, 3, 1> >, boost::range_mutable_iterator<const Eigen::Matrix<double, 3, 1> > >’:
/usr/include/boost/range/iterator.hpp:63:63:   required from ‘struct boost::range_iterator<const Eigen::Matrix<double, 3, 1> >’
/usr/include/boost/range/begin.hpp:112:61:   required by substitution of ‘template<class T> typename boost::range_iterator<const T>::type boost::range_adl_barrier::begin(const T&) [with T = Eigen::Matrix<double, 3, 1>]’
/usr/include/boost/numeric/odeint/algebra/range_algebra.hpp:52:45:   required from ‘static void boost::numeric::odeint::range_algebra::for_each3(S1&, S2&, S3&, Op) [with S1 = Eigen::Matrix<double, 3, 1>; S2 = const Eigen::Matrix<double, 3, 1>; S3 = const Eigen::Matrix<double, 3, 1>; Op = boost::numeric::odeint::default_operations::scale_sum2<double, double>]’
/usr/include/boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp:128:9:   required from ‘void boost::numeric::odeint::runge_kutta_dopri5<State, Value, Deriv, Time, Algebra, Operations, Resizer>::do_step_impl(System, const StateIn&, const DerivIn&, boost::numeric::odeint::runge_kutta_dopri5<State, Value, Deriv, Time, Algebra, Operations, Resizer>::time_type, StateOut&, DerivOut&, boost::numeric::odeint::runge_kutta_dopri5<State, Value, Deriv, Time, Algebra, Operations, Resizer>::time_type) [with System = main()::__lambda0; StateIn = Eigen::Matrix<double, 3, 1>; DerivIn = Eigen::Matrix<double, 3, 1>; StateOut = Eigen::Matrix<double, 3, 1>; DerivOut = Eigen::Matrix<double, 3, 1>; State = Eigen::Matrix<double, 3, 1>; Value = double; Deriv = Eigen::Matrix<double, 3, 1>; Time = double; Algebra = boost::numeric::odeint::range_algebra; Operations = boost::numeric::odeint::default_operations; Resizer = boost::numeric::odeint::initially_resizer; boost::numeric::odeint::runge_kutta_dopri5<State, Value, Deriv, Time, Algebra, Operations, Resizer>::time_type = double]’
/usr/include/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp:167:9:   required from ‘typename boost::disable_if<boost::is_same<StateInOut, Time>, void>::type boost::numeric::odeint::explicit_error_stepper_fsal_base<Stepper, Order, StepperOrder, ErrorOrder, State, Value, Deriv, Time, Algebra, Operations, Resizer>::do_step(System, const StateIn&, boost::numeric::odeint::explicit_error_stepper_fsal_base<Stepper, Order, StepperOrder, ErrorOrder, State, Value, Deriv, Time, Algebra, Operations, Resizer>::time_type, StateOut&, boost::numeric::odeint::explicit_error_stepper_fsal_base<Stepper, Order, StepperOrder, ErrorOrder, State, Value, Deriv, Time, Algebra, Operations, Resizer>::time_type) [with System = main()::__lambda0; StateIn = Eigen::Matrix<double, 3, 1>; StateOut = Eigen::Matrix<double, 3, 1>; Stepper = boost::numeric::odeint::runge_kutta_dopri5<Eigen::Matrix<double, 3, 1> >; short unsigned int Order = 5u; short unsigned int StepperOrder = 5u; short unsigned int ErrorOrder = 4u; State = Eigen::Matrix<double, 3, 1>; Value = double; Deriv = Eigen::Matrix<double, 3, 1>; Time = double; Algebra = boost::numeric::odeint::range_algebra; Operations = boost::numeric::odeint::default_operations; Resizer = boost::numeric::odeint::initially_resizer; typename boost::disable_if<boost::is_same<StateInOut, Time>, void>::type = void; boost::numeric::odeint::explicit_error_stepper_fsal_base<Stepper, Order, StepperOrder, ErrorOrder, State, Value, Deriv, Time, Algebra, Operations, Resizer>::time_type = double]’
/home/tlycken/exjobb/Code/alpha-orbit-follower/test/algebra/algebra-tests.cpp:21:137:   required from here
/usr/include/boost/mpl/eval_if.hpp:60:31: error: no type named ‘type’ in ‘boost::mpl::eval_if_c<true, boost::range_const_iterator<Eigen::Matrix<double, 3, 1> >, boost::range_mutable_iterator<const Eigen::Matrix<double, 3, 1> > >::f_ {aka struct boost::range_const_iterator<Eigen::Matrix<double, 3, 1> >}’
     typedef typename f_::type type;

I tried to dig into the Boost header where the error occurs, but I didn't understand enough of what's going on to be able to fix my code. Since the odeint library documentation clearly states

The main focus of odeint is to provide numerical methods implemented in a way where the algorithm is completely independent on the data structure used to represent the state x.

I believe this shouldn't be too hard to get working even if odeint doesn't support Eigen natively.

Latitudinarian

And before anyone asks: I prefer Eigen's matrix class over the built-in ones in Boost or std because of its linear algebra capabilities. I need to treat my state vector as an actual vector in state-space, where operations like vector subtraction, norms and dot products are supported. I hoped using Eigen::Matrix would be simpler than rolling my own wrapper around std::array... – Latitudinarian

Can you try to replace the definition of the stepper to runge_kutta_dopri5<state,double,state,double,vector_space_algebra> stepper; And I think you also need to include the vector space algebra: #include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>. – Anglofrench

@headmyshoulder: Thanks! I just found that out from this example and was trying it out when I got back here and saw your comment. If I make that change, my code does indeed compile, although the point it prints is the original starting point rather than the next one. I think that might be an error somewhere else (between the screen and keyboard) though... Will get back after more troubleshooting! – Latitudinarian

You need to pass dxdt as reference. Odeint expects that the r.h.s. of the ODE is stored in this variable, therefore, your need to pass it as reference: [](const state &x, state &dxdt, const double t) -> void { dxdt = x; } – Anglofrench

Ha - error between keyboard and screen, then =) If you'd like to type up an answer, I can give you rep for the help. Otherwise I'll do it myself in a while, for the next one stumbling here =) – Latitudinarian

I will do it in a second :) – Anglofrench

Answer 1

You only need to replace the definition of the stepper by

runge_kutta_dopri5<state,double,state,double,vector_space_algebra> stepper;

Eigen should work out of the Box with the vector_space_algebra but you need to specify them manually. In the next odeint version we have a mechanism for automatically detecting the algebra.

Btw. youe definition of the ODE is not correct, you need a reference for the derivatve

stepper.do_step([](const state &x, state &dxdt, const double t) -> void { 
    dxdt = x;
}, X0, 0.0, xout, 0.01);

Anglofrench

Thanks a lot! Guessing from your username, I think you might be the right person to tell about this ;) I did figure this out on my own, eventually, thanks to a small comment at the bottom of an example file in an example that was not listed in any TOC, but only in the "all examples" list (this one). With a feature that is fronted to this extent in the documentation (named as the focus of the library on the overview page) it would be nice if the documentation on how to achieve this goal was easier to reach =) – Latitudinarian

Answer 2

This solution does not seem to work for adaptive integrator

typedef runge_kutta_dopri5<state,double,state,double,vector_space_algebra> stepper_type;
auto rhs = [](const state &x, state &dxdt, const double t) -> void { dxdt = x;}
integrate_adaptive( make_controlled<stepper_type>( 1E-12 , 1E-12 ), rhs , X0 , 0. , xout , 0.01;)

Reeducate

Have a look at the boost odeint doc for adaptive integration, the state type needs to support element wise division, and element wise abs, this is not supported for Eigen::Matrix (only for Eigen::Array) hence, you somehow need to tell boost how to do that with Eigen::Matrix (transform to array and backwards)! – Stapes

I had a similar problem using Eigen::VectorXd types for state. The solution was to use #include <boost/numeric/odeint/external/eigen/eigen_algebra.hpp> (got that from here) – Calamite