Собственные векторы как ODEINT интегрируют параметры - PullRequest
Новая
       44

Собственные векторы как ODEINT интегрируют параметры

0 голосов
/ 14 апреля 2020

Я пытаюсь перенести учебник Осциллятора Harmoni c из ODEINT в Eigen, чтобы я мог использовать VectorXd для векторов состояния. Однако я не могу заставить его скомпилироваться.

Я следил за некоторыми вопросами , например , это самое близкое за исключением того, что я не используйте здесь степпер.

Это код: 

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

Компилятор жалуется на No matching function for call to 'begin' в range_algebra.hpp из ODEINT в классах и функциях integrate. На самом деле матрицы Eigen не имеют begin / end.

Я пытался поиграть с параметрами шаблона (как вы видите) безрезультатно.

Any подсказка?

Утверждение не удалось с помощью Eigen из репо

Используя последний 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
    }
  }
}
...