【问题标题】：Eigen Vectors as ODEINT integrate parameters特征向量作为 ODEINT 积分参数
【发布时间】：2020-07-27 10:38:57
【问题描述】：

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

我一直在关注some questions，例如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;
}

编译器在 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) $1 = {
  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) $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
    }
  }
}

【问题讨论】：

试试Eigen的master分支，begin()/end()前段时间已经介绍过了。
对boost::odeintStateType的要求在here中有详细说明； TL;DR：对于非固定大小的类型，除了提供begin() 和end() 之外，您的类型还必须支持boost::size( x ) 和x.resize( boost::size( y ) )
自上次发布以来已经有一段时间了——这主要是维护者的时间问题。您当然可以编写一个包装器，它只返回 this->data() 为 .begin() 和 this->data()+this->size() 为 end() （加上 const/non-const 变体和一些 typedef）
如果这在您的代码中失败，您可以只使用dxdt.resize(2)（不知道为什么它可以与std::vector 一起使用，但在这种情况下不能与Eigen::VectorXd 一起使用。但实际上，如果大小已知在编译时为 2，只需使用 Eigen::Vector2d。
您在设置成员之前尝试过dxdt.resize(2) 或dxdt.resizeLike(x) 吗？

标签： c++ boost eigen odeint

【解决方案1】：

我发现 ODEINT 实际上可以很好地与 Eigen 配合使用...只是它没有记录，据我所知。

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

你瞧，它完美无瑕：

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

希望对他人有所帮助。

【讨论】：