作为 ODEINT 积分参数的特征向量

Question

我正在尝试将谐波振荡器教程从 ODEINT 移植到 Eigen，以便我可以使用 VectorXd 作为状态向量。但是，我无法编译它。

我一直在关注 some questions, for instance this is the closest 除了我在这里不使用任何步进器。

这是代码：

#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;
}

编译器在 class 和 integrate 的函数版本中抱怨来自 ODEINT 的 range_algebra.hpp 中的 No matching function for call to 'begin'。事实上，特征矩阵没有begin/end.

我试过使用模板参数（如您所见）但无济于事。

有什么提示吗？

使用 repo 中的 Eigen 断言失败

使用 repo 中的最新 Eigen（不是最新的稳定版本），我可以按照建议编译它并运行。但是，它 使 integrate 调用树中的断言 失败：

Assertion failed: (index >= 0 && index < size()), function operator(), file eigen/include/eigen3/Eigen/src/Core/DenseCoeffsBase.h, line 427.

失败的调用是 dxdt(0) = x(1);，随后在 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);
    }

ODEINT 是否可能试图默认构造一个 VectorXd？我按照我的 ODE 调用路径进行操作， dxdt 实际上是 NULL:

(lldb) e dxdt
(state_type)  = {
  Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > = {
    m_storage = {
      m_data = 0x0000000000000000
      m_rows = 0
    }
  }
}

更糟糕的是，当使用 resizeLike 允许调整 dxdt 大小时，在第二步（因此 integrate 的第一个真正计算）我有一个 x NULL 个值：

(lldb) e dxdt
(state_type) [=14=] = {
  Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > = {
    m_storage = {
      m_data = 0x0000000000000000
      m_rows = 0
    }
  }
}
(lldb) e x
(state_type)  = {
  Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > = {
    m_storage = {
      m_data = 0x0000000000000000
      m_rows = 0
    }
  }
}

Answer 1

我发现 ODEINT 实际上可以很好地与 Eigen 一起工作...只是据我所知没有记录。

深入挖掘 ODEINT 的代码，我在 external 目录的深处发现了一个很有前途的 eigen.hpp header。

你瞧，它完美地工作：

#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;
}

希望对其他人有所帮助。

作为 ODEINT 积分参数的特征向量

Eigen Vectors as ODEINT integrate parameters

c++

boost

eigen

odeint

使用 repo 中的 Eigen 断言失败