Eigen::VectorXd and Boost::Odeint, not working

Question

I am testing with Eigen::VectorXd as state_type for boost::odeint. I use this code:

#include <Eigen/Eigen>
#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/external/eigen/eigen_algebra.hpp>

#include <iostream>
#include <vector>

template<class T>
struct push_back_state_and_time
{
    std::vector<T>& m_states;
    std::vector< double >& m_times;

    push_back_state_and_time( std::vector<T> &states ,std::vector<double> &times )
    : m_states(states) , m_times(times) { }

    void operator()(const T &x ,double t )
    {
        m_states.push_back(x);
        m_times.push_back(t);
    }
};

template<class T>
struct write_state
{
    void operator() ( const T &x, const double t ) {
        std::cout << t << "\t";
        for(size_t i = 0; i < x.size(); ++i)
            std::cout << x[i] << "\t";
        std::cout << std::endl;
    };
};

template<class T>
class odeClass {
private:
    double Ka, Kel, Vd;
public:
    odeClass(double ka, double kel, double vd) : Ka(ka), Kel(kel), Vd(vd) {};
    void operator() (const T &x, T &dxdt, const double t) {
        dxdt[0] = - Ka * x[0];
        dxdt[1] = Ka * x[0] - Kel * x[1];
    };
};

void testODE_Eigen() {
    double Ka = 0.195, Vd = 13.8, Kel = 0.79 / Vd;
    Eigen::VectorXd x(2);
    x << 40 / Vd, 0.0;

    odeClass<Eigen::VectorXd> myClass(Ka, Kel, Vd);

    boost::numeric::odeint::runge_kutta4<Eigen::VectorXd, double, Eigen::VectorXd, double, boost::numeric::odeint::vector_space_algebra> stepper;

    size_t steps = integrate_adaptive( stepper, myClass, x ,0.0 ,100.0 ,0.5 ,write_state<Eigen::VectorXd>() );
}

void testODE_Vector() {
    double Ka = 0.195, Vd = 13.8, Kel = 0.79 / Vd;
    std::vector<double> x = { 40 / Vd, 0.0 };

    odeClass<std::vector<double>> myClass(Ka, Kel, Vd);

    boost::numeric::odeint::runge_kutta4<std::vector<double>> stepper;

    size_t steps = integrate_adaptive( stepper, myClass, x ,0.0 ,100.0 ,0.5 ,write_state<std::vector<double>>() );
}

int main()
{
    testODE_Eigen();
    return 0;
}

When running the function testODE_Vector();, it works perfectly, but when runningtestODE_Eigen();` I get the first integration step, one assertion stop: "Assertion failed: index >= 0 && index < size(), file C:\Eigen-3.3.7\Eigen\src\Core\DenseCoeffsBase.h, line 408.", and a windows message saying that the application has stop working, with no clue, if a use Code::Blocks. If I run the same on Visual Studio 2017 console application, I get one error saying Cannot write on a memory address without printing anything.

Any clues?

Thank you.

Odeint does not support Eigen containers as state types. It used to work, but some upgrade - either Eigen or boost or both, broke the compatibilty a few years ago. The best solution is to use std::vector. You can use Eigen::Map to efficiently switch between Egen::VectorXd and std::vector. — RHertel

torrey1028 torrey1028 · Accepted Answer · 2020-06-05T22:58:14

It looks like you are failing an assertion inside the Eigen library you are using. With a message like Assertion failed: index >= 0 && index < size() the library is probably trying to iterate over a vector internally and checks that the vector is valid before trying to access it. I would check the objects you are passing in and make sure they are valid.

It looks like one of the main differences in the two function testODE_Vector() and testODE_Eigen() is the way that you create that x. I'm not sure what this line x << 40 / Vd, 0.0; is intended to do, but I would start there and verify that the value of x is right before it's passed into integrate_adaptive

Eigen::VectorXd and Boost::Odeint, not working

1 Answers