While trying to mimic a Python Code in CPP (With Eigen)

Python Code :

import numpy as np 
from numpy.linalg import inv
from matplotlib import pyplot as plt
from math import cos, sin, pi


def G(y,t): 
    a1d, a2d = y[0], y[1]
    a1, a2 = y[2], y[3]

    m11, m12 = (m1 m2)*l1, m2*l2*cos(a1-a2)
    m21, m22 = l1*cos(a1-a2), l2
    m = np.array([[m11, m12],[m21, m22]])

    f1 = -m2*l2*a2d*a2d*sin(a1-a2) - (m1 m2)*g*sin(a1)
    f2 = l1*a1d*a1d*sin(a1-a2) - g*sin(a2)
    f = np.array([f1, f2])

    accel = inv(m).dot(f)

    return np.array([accel[0], accel[1], a1d, a2d])


def RK4_step(y, t, dt):
    k1 = G(y,t)
    k2 = G(y 0.5*k1*dt, t 0.5*dt)
    k3 = G(y 0.5*k2*dt, t 0.5*dt)
    k4 = G(y k3*dt, t dt)

    return dt * (k1   2*k2   2*k3   k4) /6

# variables
m1, m2 = 2.0, 1.0
l1, l2 = 1.0, 2.0
g = 9.81

delta_t = 0.01
time = np.arange(0.0, 10.0, delta_t)

# initial state
y = np.array([0,0,0,1.0])   # [velocity, displacement]

Y1 = []
Y2 = []

# time-stepping solution
for t in time:
    y = y   RK4_step(y, t, delta_t) 

    Y1.append(y[2])
    Y2.append(y[3])


# plot the result
plt.plot(time,Y1)
plt.plot(time,Y2)
plt.grid(True)
plt.legend(['Angle1', 'Angle2'], loc='lower right')
plt.show()

Source : https://github.com/apf99/Double-Pendulum/blob/master/d_pen.py (https://www.youtube.com/watch?v=nBBQKZb6JZk)

Being a Noob to C and Eigen.

Transcribed Code in C :

#include<iostream>
#include<math.h>
#include<Eigen/Eigen>
#include<fstream>
// grc : gravitational constant

#define grc 9.81


// solving for acceleration
void G(const Eigen::Ref<const Eigen::VectorXd>& dummy2, const double& t, const double& m1, const double& m2, const double& l1, const double& l2, Eigen::Ref<Eigen::Vector4d>& k){
  // angular velocities (av1, av2) and displacements (a1, a2),
  double av1, av2, a1, a2;
  // angular velocities
  av1 = dummy2[0];
  av2 = dummy2[1];
  // angular position
  a1 = dummy2[3];
  a2 = dummy2[4];

  // Defiing the mass like matrix
  Eigen::Matrix2d A;
  A<<(m1 m2)*l1, m2*l2*std::cos(a1-a2),  l1*std::cos(a1-a2), l2;
  // Defining the forcing like function
  Eigen::Vector2d F, acceleration;
  F<<-m2 * l2 * av2*av2*std::sin(a1-a2) - (m1 m2)*grc*std::sin(a1), l1*av1*av1*std::sin(a1-a2) - grc*std::sin(a2);

  // product of the inverse mass and forcing function
  acceleration = (A.inverse()) * F;

  //assiging the values
  k(0) = acceleration(0);
  k(1) = acceleration(1);
  k(2) = av1;
  k(3) = av2;
   
}

// attempt to pass by refernce
// Ranga Kutta 4
void RK4(const Eigen::Ref<Eigen::VectorXd>& dummy, const double& t, const double& dt, const double& m1, const double& m2, const double& l1, const double& l2){
  Eigen::Vector4d k1 , k2 , k3 , k4 ;
  G(dummy, t, m1, m2, l1, l2, k1);
  G(dummy  0.5 * dt * k1, t 0.5*dt, m1, m2, l1, l2, k2);
  G(dummy  0.5 * dt * k2, t 0.5*dt, m1, m2, l1, l2, k3);
  G(dummy  dt * k3, t dt, m1, m2, l1, l2, k4);

  dummy = (dt/6.0) * (k1   2*k2   2* k3   k4);
}


int main(int argc, char** argv){

  // m1 and m2 are the masses of the pendulum 
  double m1 = 2.0, m2 = 1.0;
  // l1 and l2 are the lengths of the masses
  double l1 = 1.0, l2 = 2.0;
  // time delta
  double dt = 0.01;
  // Defining the angular displacement y(2), y(3) and  angular velocities y(0), y(1)
  Eigen::VectorXd y(4), y_dummy(4);
  y(0) = 0.0;
  y(1) = 0.0;
  y(2) = 0.0;
  y(3) = 1.0;

  // time as t
  double t = 0.0;

  //file opening to write data
  std::fstream pend("double_pendulum.dat",std::ios::out|std::ios::app);
  //position and velocity are written to the dat file
  pend<<"v(x)  v(y) p(x) p(y) \n";
  
  // time ; stepping
  for(unsigned i = 0; i<1000;i  ){
    t  = i/100;
    y_dummy = y;
    RK4(y_dummy,t,dt, m1, m2, l1, l2);
    y  = y_dummy;
    pend<<(y.transpose())<<std::endl;
  }

  pend.close();
  return 0;
  
}

When compiled with g -std=c 11 -I /usr/include/eigen3 -o dp dp.cpp

The following errors popup (The jist of the error is mentioned below)

error: cannot bind non-const lvalue reference of type ‘Eigen::Ref<Eigen::Matrix<double, 4, 1> >&’ to an rvalue of type ‘Eigen::Ref<Eigen::Matrix<double, 4, 1> >’
   43 |   G(dummy, t, m1, m2, l1, l2, k1);

etc .. etc.. etc.

for each non const reference and matrix passed.

When either or both are defined as const, they become immutable, hence the vector will not be updated.

Being completely flabbergasted by Eigen_tux documentation , what could be the solution in layman's terms.

Or, are there libraries that can perform this task in a simple manner?

PS : Couldn't compile the script completely

CodePudding user response：

A mutable reference, more precisely an l-value reference, only ever points at an existing, named object (the actual language rules are of course more complicated but that is the intention). The Eigen::Ref object is not one of those since it is created as a temporary object in the function call. Therefore the compiler refuses to create such a reference.

As a consequence, it is customary to pass inputs as const Eigen::Ref<const Matrix>& and input/outputs as Eigen::Ref<Matrix>. The function signature should look like this:

void G(const Eigen::Ref<const Eigen::VectorXd>& dummy2,
    double t, double m1, double m2, double l1, double l2,
    Eigen::Ref<Eigen::Vector4d> k);

Note that the Ref object functionally references the original matrix or matrix block. It doesn't need to be a reference in the C sense to do this. You could also pass the input Eigen::Ref<const Matrix> by value. However, this would introduce a slight loss in performance if you continue to pass that Ref to a sub-function, since it would create a new copy.

For mutable references we accept that cost for consistency and because mutable Refs are actually cheaper objects than Ref<const Matrix>. The reason being that const-Refs copy their inputs if the inner stride would not be 1, e.g. when you pass a single row of a matrix as a Ref<const Vector>. This is described in the Ref documentation.

Note also that passing simple doubles as const reference is pointless. Pass-by-value for scalars is more efficient.