使用Boost :: odeint與特徵::矩陣作爲狀態向量

Question

我試圖利用ODE integration capabilities of Boost使用Matrix class from Eigen 3作爲我的狀態向量，但我遇到了深入到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; 
} 

如果我COMENT了呼叫stepper.do_step一切編譯和運行得很好，但當然沒有做什麼有趣的事。如果我不這樣做，升壓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，我相信這應該不會太難。

Answer 1

你只需要通過

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

徵應開箱與vector_space_algebra更換步進的定義，但需要手動指定它們。在接下來的版本中，我們有一個自動檢測代數的機制。

Btw。悠逸的ODE的定義是不正確的，你需要的derivatve

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

Answer 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;)