Использование std :: vector <Eigen :: Vector3d> как тип состояния в odeint

Question

В настоящее время я пытаюсь использовать odeint и Eigen3 для интеграции системы nBody (цель - это библиотека, предоставляющая расширенные подпрограммы для использования в образовании планет, такие как как смешанная переменная симплектическая или вариационная вариация МВМС МВС). Экспериментируя с различными степпер я обнаружил, что при использовании нормальных степпер 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?

и второй:

/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, так как наличие сгруппированных координат очень полезно для удобочитаемости правых сторон ОДУ.

Вот приведенный пример кода, который я использовал.

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

Головной уборный раствор работал. Я создал алгебру, унаследованный от 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, но я еще не пробовал.

Answer 1

Я боюсь, что ваш тип состояния может быть легко интегрирован. Проблема состоит в том, что смешивается тип типа типа (vector<>) с типами состояний, подобными вектору (Vector3d). Вам понадобится range_algebra, чтобы перебрать вектор. Но norm_inf от range_algebra не работает в вашем случае.

Что вы можете сделать, это дублировать range_algebra и оставить все методы for_eachX такими, какие они есть, и только повторно реализовать norm_inf в соответствии с вашим типом состояния. Это не должно быть сложно. Затем вы используете свою новую алгебру просто через

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