Search code examples
c++matlabeigeneigen3odeint

Using Eigen with Odeint to do matrix*vector multiplication inside dxdt function


I'm trying to write some code in C++ that would otherwise be quite easy to write in MatLab, using Eigen and Odeint. However, I am new to both the Eigen and Odeint libraries and I am not getting very far. If someone out there could just point me in the right direction (in terms of how to code it), I should be able to take it from there.

Here's roughly what I'd like to do, but can't get to work (in pseudo code; anything in all caps is a constant):

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

Many thanks to anyone who can lend a hand.


Solution

  • The (Eigen) syntax problems with your code are as follows:

    1. v.colwise.sum() should be v.colwise().sum()
    2. Likewise, v.rowwise.sum() should be v.rowwise().sum()
    3. state_type_1d u(N,N); is mixed vector type, matrix constructor.
    4. v = (0.5 * (1 + tanh((x-0.5) * GAIN))); appears to be a coefficient wise vector product into a matrix variable. And I don't think that Eigen has a built in tanh but I could be wrong. Even if so, wrong syntax. It should be (x-0.5).tanh() if it exists.