在 odeint 中使用 std::vector<Eigen::Vector3d> 作为状态类型答案

【问题标题】：Using a std::vector<Eigen::Vector3d> as state type in odeint在 odeint 中使用 std::vector<Eigen::Vector3d> 作为状态类型
【发布时间】：2015-05-12 00:10:47
【问题描述】：

我目前正在尝试使用 odeint 和 Eigen3 来集成 nBody 系统（目标是提供用于行星形成的高级例程的库，例如 MVS 的混合变量辛或钱伯斯混合变体）。在尝试使用不同的步进器时，我发现使用普通步进器时 state_type std::vector<Eigen::Vector3d> 工作正常，但使用受控步进器（例如 burlisch_stoer）编译失败，第一条错误消息是：

/usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:87:40: error: cannot convert ‘boost::numeric::odeint::norm_result_type<std::vector<Eigen::Matrix<double, 3, 1> >, void>::type {aka Eigen::Matrix<double, 3, 1>}’ to ‘boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>::value_type {aka double}’ in return
    return algebra.norm_inf( x_err );

这是否意味着 norm_result_type 推导不正确？这个规范到底是做什么的？应该是 x_err 中找到的最高 value_type 吗？

第二个：

/usr/include/boost/numeric/odeint/algebra/range_algebra.hpp:132:89: error: call of overloaded ‘Matrix(int)’ is ambiguous
    static_cast< typename norm_result_type<S>::type >( 0 ) );

我必须提供自己的代数才能以这种方式使用它吗？我宁愿不切换到 std::vector<double> 或 Eigen::VectorNd，因为将坐标分组对 ODE 右侧的可读性非常有益。

这是我使用的代码的简化示例。

#include "boost/numeric/odeint.hpp"
#include "boost/numeric/odeint/external/eigen/eigen.hpp"
#include "Eigen/Core"

using namespace boost::numeric::odeint;

typedef std::vector<Eigen::Vector3d> state_type;

struct f
{
    void operator ()(const state_type& state, state_type& change, const  double /*time*/)
    {
    };
};

int main()
{
    // Using this compiles
    typedef euler <state_type> stepper_euler;
    // Using this does not compile
    typedef bulirsch_stoer <state_type> stepper_burlisch;

    state_type x;
    integrate_const(
        stepper_burlisch(),
        f(),
        x,
        0.0,
        1.0,
        0.1
        );

    return 0;
}

Headmyshoulders 解决方案有效。我创建了一个继承自range_algebra 的代数：

class custom_algebra : public boost::numeric::odeint::range_algebra {
    public:
        template< typename S >
        static double norm_inf( const S &s )
        {
            double norm = 0;
            for (const auto& inner : s){
                const double tmp = inner.maxCoeff();
                if (tmp > norm)
                    norm = tmp;
            }
            return norm;
        }
} ;

我认为应该可以使用range_algebra 和vector_space_algebra 通用地创建这样一个代数，但我还没有尝试过。

【问题讨论】：

标签： c++ vector eigen3 odeint

【解决方案1】：

我担心您的状态类型不容易被集成。问题是，它混合了类范围状态类型（vector<>）和类向量空间状态类型（Vector3d）。您需要 range_algebra 来迭代向量。但是来自range_algebra 的norm_inf 不适用于您的情况。

您可以做的是复制range_algebra 并保留所有for_eachX 方法，只重新实现norm_inf 以适应您的状态类型。那应该不难。然后你只需通过

来使用你的新代数

typedef bulirsch_stoer< state_type , double , state_type , double , your_algebra > stepper_burlisch

【讨论】：