Using a std::vector<Eigen::Vector3d> as state type in odeint

Question

I'm currently trying to use odeint and Eigen3 to integrate an nBody system (target is a library providing advanced routines for use in planet formation, such as mixed variable symplectic or Chambers hybrid variation of MVS). While experimenting with different steppers I found that when using normal steppers the state_type std::vector<Eigen::Vector3d> works fine, but with controlled steppers (such as burlisch_stoer) compilation fails, the first error message being:

/usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:87:40: error: cannot convert ‘boost::numeric::odeint::norm_result_type<std::vector<Eigen::Matrix<double, 3, 1> >, void>::type {aka Eigen::Matrix<double, 3, 1>}’ to ‘boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>::value_type {aka double}’ in return
    return algebra.norm_inf( x_err );

Does this mean that the norm_result_type is deduced incorrectly? What does this norm do exactly? Should it be the highest value_type found in x_err?

and the second:

/usr/include/boost/numeric/odeint/algebra/range_algebra.hpp:132:89: error: call of overloaded ‘Matrix(int)’ is ambiguous
    static_cast< typename norm_result_type<S>::type >( 0 ) );

Do I have to provide my own algebra to use it this way? I would rather not switch to a std::vector<double> or an Eigen::VectorNd, as having the coordinates grouped is very beneficial to the readability of the right hand sides of the ODEs.

Here is the a reduced example of the code I used.

#include "boost/numeric/odeint.hpp"
#include "boost/numeric/odeint/external/eigen/eigen.hpp"
#include "Eigen/Core"

using namespace boost::numeric::odeint;

typedef std::vector<Eigen::Vector3d> state_type;

struct f
{
    void operator ()(const state_type& state, state_type& change, const  double /*time*/)
    {
    };
};

int main()
{
    // Using this compiles
    typedef euler <state_type> stepper_euler;
    // Using this does not compile
    typedef bulirsch_stoer <state_type> stepper_burlisch;

    state_type x;
    integrate_const(
        stepper_burlisch(),
        f(),
        x,
        0.0,
        1.0,
        0.1
        );

    return 0;
}

Headmyshoulders solution worked. I have created a algebra inheriting from the range_algebra:

class custom_algebra : public boost::numeric::odeint::range_algebra {
    public:
        template< typename S >
        static double norm_inf( const S &s )
        {
            double norm = 0;
            for (const auto& inner : s){
                const double tmp = inner.maxCoeff();
                if (tmp > norm)
                    norm = tmp;
            }
            return norm;
        }
} ;

I think it should be possible to create such an algebra generically using both range_algebra and vector_space_algebra, but i haven't tried yet.

Answer 1

I fear that your state type can is not easily be integrated. The problem is, that is mixes a range-like state type (the vector<>) with vector-space-like state types (the Vector3d). You need the range_algebra to iterate over the vector. But norm_inf from range_algebra does not work in your case.

What you can do is to duplicate the range_algebra and and leave all for_eachX methods as they are and only re-implement norm_inf to fit with your state type. That should not be to difficult. Then you use your new algebra simply via

typedef bulirsch_stoer< state_type , double , state_type , double , your_algebra > stepper_burlisch