I am trying to port the Harmonic Oscillator tutorial from ODEINT to Eigen, so that I could use VectorXd for state vectors. I cannot, however, make it compile.
I've been following some questions, for instance this is the closest except that I don't use any stepper here.
This is the code:
#include <iostream>
#include <vector>
#include <Eigen/Dense>
#include <boost/numeric/odeint.hpp>
typedef Eigen::VectorXd state_type;
// was vector<double>
const double gam = 0.15;
/* The rhs of x' = f(x) defined as a class */
class harm_osc
{
public:
void operator() ( const state_type &x , state_type &dxdt , const double /* t */ )
{
dxdt(0) = x(1);
dxdt(1) = -x(0) - gam*x(1);
// dxdt[0] = x[1];
// dxdt[1] = -x[0] - gam*x[1];
}
};
/* The rhs of x' = f(x) */
void harmonic_oscillator(const state_type &x, state_type &dxdt, const double /* t */ )
{
dxdt(0) = x(1);
dxdt(1) = -x(0) - gam*x(1);
// dxdt[0] = x[1];
// dxdt[1] = -x[0] - gam*x[1];
}
void printer(const state_type &x , const double t)
{
// std::cout << t << "," << x[0] << "," << x[1] << std::endl;
std::cout << t << "," << x(0) << "," << x(1) << std::endl;
};
int main(int argc, const char * argv[])
{
state_type x(2);
x(0) = 1.0;
x(1) = 0.0;
// x[0] = 1.0;
// x[1] = 0.0;
std::cout << ">>>> FUNCTION" << std::endl;
// boost::numeric::odeint::integrate(harmonic_oscillator, x, 0.0, 10.0, 0.1, printer);
// boost::numeric::odeint::runge_kutta4<state_type, double, state_type, double, boost::numeric::odeint::vector_space_algebra> stepper();
boost::numeric::odeint::integrate<double, decltype(harmonic_oscillator), state_type, double, decltype(printer)>(harmonic_oscillator, x, 0.0, 10.0, 0.1, printer);
std::cout << ">>>> CLASS" << std::endl;
x(0) = 1.0;
x(1) = 0.0;
// x[0] = 1.0;
// x[1] = 0.0;
harm_osc ho;
boost::numeric::odeint::integrate<double, decltype(harmonic_oscillator), state_type, double, decltype(printer)>(ho, x, 0.0, 10.0, 0.1, printer);
return 0;
}
The compiler complains about No matching function for call to 'begin' in range_algebra.hpp from ODEINT in both class and function versions of integrate. As a matter of fact, Eigen matrices have no begin/end.
I've tried to play with the template parameters (as you see) with no avail.
Any hint?
Using the latest Eigen from the repo (not the latest stable version), I can, as suggested, compile it and run. However, it fails an assertion in the integrate call tree:
Assertion failed: (index >= 0 && index < size()), function operator(), file eigen/include/eigen3/Eigen/src/Core/DenseCoeffsBase.h, line 427.
The call that fails is dxdt(0) = x(1); and subsequently in DenseCoeffsBase.h:
EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE Scalar&
operator()(Index index)
{
eigen_assert(index >= 0 && index < size()); // <---- INDEX IS 0, SIZE IS 0
return coeffRef(index);
}
Is it possible that ODEINT is trying to default-construct a VectorXd? I followed the path to my ODE call and dxdt is actually NULL:
(lldb) e dxdt
(state_type) $1 = {
Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > = {
m_storage = {
m_data = 0x0000000000000000
m_rows = 0
}
}
}
What is worse is that when using resizeLike to allow resizing dxdt, in the second step (so the first real computation of integrate) I have a x with NULL values:
(lldb) e dxdt
(state_type) $0 = {
Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > = {
m_storage = {
m_data = 0x0000000000000000
m_rows = 0
}
}
}
(lldb) e x
(state_type) $1 = {
Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > = {
m_storage = {
m_data = 0x0000000000000000
m_rows = 0
}
}
}
I found that ODEINT actually works fine with Eigen... only it is not documented, as far as I can see.
Digging around ODEINT's code, I have found a promising eigen.hpp header deep in the external directory.
And lo and behold, it works flawlessly:
#include <fstream>
#include <iostream>
#include <vector>
#include <boost/numeric/odeint/external/eigen/eigen.hpp>
#include <Eigen/Eigenvalues>
#define FMT_HEADER_ONLY
#include "fmt/core.h"
#include "fmt/format.h"
#include "fmt/format-inl.h"
#include "fmt/printf.h"
using namespace std;
int main(int argc, char *argv[])
{
Eigen::VectorXd v;
v.resize(2);
typedef Eigen::VectorXd state_type;
const double gam = 0.15;
v(0) = 1.0;
v(1) = 1.1;
auto harmonic_oscillator = [&](const state_type &x, state_type &dxdt, const double t)
{
dxdt[0] = x[1];
dxdt[1] = -x[0] - gam*x[1];
};
auto printer = [&](const state_type &x, const double t)
{
fmt::print(out, "time: {} state: {}\n", t, x);
};
odeint::integrate(harmonic_oscillator, v, 0.0 , 10.0 , 0.01, printer);
return 0;
}
Hope it helps others.