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.

Answer 1

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

Answer 2

My answer is a little late but in case someone else runs into this issue, here's what I found.

The issue seems to be that OdeInt can't handle properly the dynamic size with Eigen vectors and matrices. Therefore when creating dxdt, it creates an empty dynamic matrix or vector. This leads to an issue in your operator overload, where you try to access an element of dxdt where it contains none.

A quick fix I found was to use the resize() function (or conservativeResize()) to make sure dxdt has the proper size:

void operator() (const T &x, T &dxdt, const double t) {
    dxdt.resize(x.size())
    dxdt[0] = - Ka * x[0];
    dxdt[1] = Ka * x[0] - Kel * x[1];
};

Note that if you want to use matrices instead of vectors you will have to use x.rows() and x.cols() instead of x.size().