Использование Eigen с Odeint для умножения матрицы на вектор внутри функции dxdt

Я пытаюсь написать некоторый код на С++, который в противном случае было бы довольно легко написать в MatLab, используя Eigen и Odeint. Однако я новичок в библиотеках Eigen и Odeint, и я не очень далеко продвинулся. Если бы кто-то там мог просто указать мне правильное направление (с точки зрения того, как его кодировать), я мог бы взять его оттуда.

Вот примерно то, что я хотел бы сделать, но не могу приступить к работе (в псевдокоде; все, что написано заглавными буквами, является константой):

typedef Eigen::Matrix<double, Dynamic, 1> state_type_1d;
typedef Eigen::Matrix<double, Dynamic, Dynamic> state_type_2d;

state_type Q(N*N,N*N) = ... // initialize Q

void dxdt_fun( state_type_1d &x , state_type_1d &dxdt , double t )
{

    static state_type_2d v(N,N);

    v = (0.5 * (1 + tanh((x-0.5) * GAIN)));
    dxdt = -x*LAMBDA + v.colwise.sum() + v.rowwise.sum(); // needs bsxfun??
    dxdt = dxdt + Q * x; // matrix times vector
}

void main(int argc, char **argv)
{
    state_type_1d u(N,N); 
    srand((unsigned int) time(0));
    u.setRandom(); // picks random numbers from -1 to 1
    runge_kutta_dopri5<state_type_1d,double,state_type_1d,double,vector_space_algebra>stepper;

    integrate_adaptive(stepper, dxdt_fun, u, 0.0, 100.0, 0.01, write_dxdt);
}

Большое спасибо всем, кто может протянуть руку.


person Aaron Schurger    schedule 23.08.2015    source источник


Ответы (2)


Синтаксические проблемы (Eigen) с вашим кодом следующие:

  1. v.colwise.sum() должно быть v.colwise().sum()
  2. Аналогично, v.rowwise.sum() должно быть v.rowwise().sum()
  3. state_type_1d u(N,N); — смешанный векторный тип, матричный конструктор.
  4. v = (0.5 * (1 + tanh((x-0.5) * GAIN))); выглядит как векторное произведение с коэффициентом на матричную переменную. И я не думаю, что у Eigen есть встроенный tanh, но я могу ошибаться. Даже если это так, неправильный синтаксис. Он должен быть (x-0.5).tanh(), если он существует.
person Avi Ginsburg    schedule 23.08.2015

Я написал решение, в котором использовал как Eigen, так и odeint. См. эту ссылку, в качестве примера используется двойной маятник.

person CroCo    schedule 02.10.2015