使用 Boost::odeint 和 Eigen::Matrix 作为状态向量答案

【问题标题】：Using Boost::odeint with Eigen::Matrix as state vector使用 Boost::odeint 和 Eigen::Matrix 作为状态向量
【发布时间】：2014-05-12 00:15:50
【问题描述】：

我正在尝试使用Matrix class from Eigen 3 作为我的状态向量来利用ODE integration capabilities of Boost，但我在 Boost 中遇到了我不知道如何解决的问题。

我正在尝试做的一个最小示例：

#include <Eigen/Core>
#include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
#include <iostream>

using namespace Eigen;
using namespace boost::numeric::odeint;

template<size_t N>
using vector = Matrix<double, N, 1>;

typedef vector<3> state;

int main() {

    state X0;
    X0 << 1., 2., 3.;
    state xout = X0;

    runge_kutta_dopri5<state> stepper;

    // If I remove these lines, everything compiles fine
    stepper.do_step([](const state x, state dxdt, const double t) -> void { 
        dxdt = x;
    }, X0, 0.0, xout, 0.01);

    std::cout << xout << std::endl;
}

如果我取消对stepper.do_step 的调用，一切都会编译并运行得很好，但当然不会做任何有趣的事情。如果我没有，Boost vomits compile errors over my terminal，第一个是

In file included from /usr/include/boost/mpl/aux_/begin_end_impl.hpp:20:0,
                 from /usr/include/boost/mpl/begin_end.hpp:18,
                 from /usr/include/boost/mpl/is_sequence.hpp:19,
                 from /usr/include/boost/fusion/support/detail/is_mpl_sequence.hpp:12,
                 from /usr/include/boost/fusion/support/tag_of.hpp:13,
                 from /usr/include/boost/fusion/support/is_sequence.hpp:11,
                 from /usr/include/boost/fusion/sequence/intrinsic_fwd.hpp:12,
                 from /usr/include/boost/fusion/sequence/intrinsic/front.hpp:10,
                 from /usr/include/boost/fusion/include/front.hpp:10,
                 from /usr/include/boost/numeric/odeint/util/is_resizeable.hpp:26,
                 from /usr/include/boost/numeric/odeint/util/state_wrapper.hpp:25,
                 from /usr/include/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp:27,
                 from /usr/include/boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp:24,
                 from /home/tlycken/exjobb/Code/alpha-orbit-follower/test/algebra/algebra-tests.cpp:2:
/usr/include/boost/mpl/eval_if.hpp: In instantiation of ‘struct boost::mpl::eval_if_c<true, boost::range_const_iterator<Eigen::Matrix<double, 3, 1> >, boost::range_mutable_iterator<const Eigen::Matrix<double, 3, 1> > >’:
/usr/include/boost/range/iterator.hpp:63:63:   required from ‘struct boost::range_iterator<const Eigen::Matrix<double, 3, 1> >’
/usr/include/boost/range/begin.hpp:112:61:   required by substitution of ‘template<class T> typename boost::range_iterator<const T>::type boost::range_adl_barrier::begin(const T&) [with T = Eigen::Matrix<double, 3, 1>]’
/usr/include/boost/numeric/odeint/algebra/range_algebra.hpp:52:45:   required from ‘static void boost::numeric::odeint::range_algebra::for_each3(S1&, S2&, S3&, Op) [with S1 = Eigen::Matrix<double, 3, 1>; S2 = const Eigen::Matrix<double, 3, 1>; S3 = const Eigen::Matrix<double, 3, 1>; Op = boost::numeric::odeint::default_operations::scale_sum2<double, double>]’
/usr/include/boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp:128:9:   required from ‘void boost::numeric::odeint::runge_kutta_dopri5<State, Value, Deriv, Time, Algebra, Operations, Resizer>::do_step_impl(System, const StateIn&, const DerivIn&, boost::numeric::odeint::runge_kutta_dopri5<State, Value, Deriv, Time, Algebra, Operations, Resizer>::time_type, StateOut&, DerivOut&, boost::numeric::odeint::runge_kutta_dopri5<State, Value, Deriv, Time, Algebra, Operations, Resizer>::time_type) [with System = main()::__lambda0; StateIn = Eigen::Matrix<double, 3, 1>; DerivIn = Eigen::Matrix<double, 3, 1>; StateOut = Eigen::Matrix<double, 3, 1>; DerivOut = Eigen::Matrix<double, 3, 1>; State = Eigen::Matrix<double, 3, 1>; Value = double; Deriv = Eigen::Matrix<double, 3, 1>; Time = double; Algebra = boost::numeric::odeint::range_algebra; Operations = boost::numeric::odeint::default_operations; Resizer = boost::numeric::odeint::initially_resizer; boost::numeric::odeint::runge_kutta_dopri5<State, Value, Deriv, Time, Algebra, Operations, Resizer>::time_type = double]’
/usr/include/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp:167:9:   required from ‘typename boost::disable_if<boost::is_same<StateInOut, Time>, void>::type boost::numeric::odeint::explicit_error_stepper_fsal_base<Stepper, Order, StepperOrder, ErrorOrder, State, Value, Deriv, Time, Algebra, Operations, Resizer>::do_step(System, const StateIn&, boost::numeric::odeint::explicit_error_stepper_fsal_base<Stepper, Order, StepperOrder, ErrorOrder, State, Value, Deriv, Time, Algebra, Operations, Resizer>::time_type, StateOut&, boost::numeric::odeint::explicit_error_stepper_fsal_base<Stepper, Order, StepperOrder, ErrorOrder, State, Value, Deriv, Time, Algebra, Operations, Resizer>::time_type) [with System = main()::__lambda0; StateIn = Eigen::Matrix<double, 3, 1>; StateOut = Eigen::Matrix<double, 3, 1>; Stepper = boost::numeric::odeint::runge_kutta_dopri5<Eigen::Matrix<double, 3, 1> >; short unsigned int Order = 5u; short unsigned int StepperOrder = 5u; short unsigned int ErrorOrder = 4u; State = Eigen::Matrix<double, 3, 1>; Value = double; Deriv = Eigen::Matrix<double, 3, 1>; Time = double; Algebra = boost::numeric::odeint::range_algebra; Operations = boost::numeric::odeint::default_operations; Resizer = boost::numeric::odeint::initially_resizer; typename boost::disable_if<boost::is_same<StateInOut, Time>, void>::type = void; boost::numeric::odeint::explicit_error_stepper_fsal_base<Stepper, Order, StepperOrder, ErrorOrder, State, Value, Deriv, Time, Algebra, Operations, Resizer>::time_type = double]’
/home/tlycken/exjobb/Code/alpha-orbit-follower/test/algebra/algebra-tests.cpp:21:137:   required from here
/usr/include/boost/mpl/eval_if.hpp:60:31: error: no type named ‘type’ in ‘boost::mpl::eval_if_c<true, boost::range_const_iterator<Eigen::Matrix<double, 3, 1> >, boost::range_mutable_iterator<const Eigen::Matrix<double, 3, 1> > >::f_ {aka struct boost::range_const_iterator<Eigen::Matrix<double, 3, 1> >}’
     typedef typename f_::type type;

我试图深入研究the Boost header where the error occurs，但我对正在发生的事情了解不够，无法修复我的代码。由于odeint库文档clearly states

odeint 的主要重点是提供以算法完全独立于用于表示状态 x 的数据结构的方式实现的数值方法。

我相信即使 odeint 本身不支持 Eigen，这也不难实现。

【问题讨论】：

在有人问之前：我更喜欢 Eigen 的矩阵类而不是 Boost 或 std 中的内置矩阵类，因为它具有线性代数功能。我需要将我的状态向量视为状态空间中的实际向量，其中支持向量减法、范数和点积等操作。我希望使用Eigen::Matrix 比在std::array 周围滚动我自己的包装更简单...
能否尝试将stepper的定义替换为runge_kutta_dopri5 stepper；而且我认为您还需要包含向量空间代数：#include .
@headmyshould：谢谢！我刚从this example 发现了这一点，当我回到这里并看到你的评论时正在尝试。如果我进行该更改，我的代码确实可以编译，尽管它打印的点是原始起点而不是下一个起点。我认为这可能是其他地方（屏幕和键盘之间）的错误......在更多故障排除后会回来！
您需要传递 dxdt 作为参考。 Odeint 预计 r.h.s. ODE 的值存储在此变量中，因此，您需要将其作为参考传递： [](const state &x, state &dxdt, const double t) -> void { dxdt = x; }
哈 - 键盘和屏幕之间的错误，然后 =) 如果您想输入答案，我可以为您提供帮助。否则我一会儿自己做，下一个在这里磕磕绊绊=）

标签： c++ templates boost eigen template-meta-programming

【解决方案1】：

只需将步进器的定义替换为

runge_kutta_dopri5<state,double,state,double,vector_space_algebra> stepper;

Eigen 应该可以使用 vector_space_algebra 开箱即用，但您需要手动指定它们。在下一个 odeint 版本中，我们有一个自动检测代数的机制。

顺便说一句。您对 ODE 的定义不正确，您需要派生的参考

stepper.do_step([](const state &x, state &dxdt, const double t) -> void { 
    dxdt = x;
}, X0, 0.0, xout, 0.01);

【讨论】：

非常感谢！从您的用户名猜测，我认为您可能是合适的人来谈论这个；）我确实自己解决了这个问题，最终，这要归功于示例文件底部的一个小注释，该示例中未列出任何 TOC，但仅在“所有示例”列表中 (this one)。有了一个在文档中如此突出的特性（在概述页面上被命名为库的焦点），如果关于如何实现这个目标的文档更容易达到，那就太好了=）

【解决方案2】：

此解决方案似乎不适用于自适应积分器

typedef runge_kutta_dopri5<state,double,state,double,vector_space_algebra> stepper_type;
auto rhs = [](const state &x, state &dxdt, const double t) -> void { dxdt = x;}
integrate_adaptive( make_controlled<stepper_type>( 1E-12 , 1E-12 ), rhs , X0 , 0. , xout , 0.01;)

【讨论】：

自适应积分看the boost odeint doc，状态类型需要支持element-wise除法和element-wise abs，Eigen::Matrix不支持这个（仅Eigen::Array ) 因此，您需要以某种方式告诉 boost 如何使用 Eigen::Matrix 执行此操作（转换为数组并向后转换）！
我在使用 Eigen::VectorXd 类型作为状态时遇到了类似的问题。解决方案是使用#include <boost/numeric/odeint/external/eigen/eigen_algebra.hpp>（从here获得）