Quick tutorial

Introduction

This tutorial is intended to quickly be able to use the basics of the library without diving into the internals. Developpers might consider reading the documentation instead.

RKMK module

In this module, we consider dynamical systems governed by the equation

\[ \dot Y(t) = A(Y)Y(t) \]

where \(Y\) lies on the configuration manifold \(\mathcal M\subset\mathbb R^n\) of the system, with \(\mathcal M\) left-invariant for the action of the Lie group \(G\), and \(A(Y)\) belongs to the associated Lie algebra \(\mathfrak g\) for all \(Y\).

To use this module, we need only need to include the RKMK header of Geomi.

#include <Geomi/RKMK>

The first step is to define the Lie algebra. This is done either by using a preexisting implementation of the base class LieAlgebraBase or by provinding your own. In this case, we use the class SO3::Algebra implementing the \(\mathfrak{so}(3)\) algebra.

typedef SO3::Algebra<double> Algebra;

The second step is to provide the equations for the problem, namely the matrix \(A\). This is done by inheriting the abstract class RKMK::Abstract::Problem and implementing the virtual functions computeA and computeJacobianA. In our example with the rigid body problem, we have

class RigidBodyProblem : public RKMK::Abstract::Problem<Algebra>
{
    private:
        // Inverted inertia tensor
        Eigen::Matrix<double,3,1> m_Iinv;
    
    public:
        /* ... */

        bool
        computeA (Algebra& A, const Q& x)
        {
            A = Algebra((-m_Iinv).asDiagonal()*x);
            return true;
        }

        bool
        computeJacobianA (std::vector<Algebra>& JA, const Q& x)
        {
            Eigen::Matrix<double,3,3> M = (-m_Iinv).asDiagonal();
            for (int i=0; i<3; i++)
                JA.push_back(Algebra(M.col(i)));
            return true;
        }
};

The system is now fully defined and can be initialized and solved using a Runge-Kutta Munthe-Kaas integrator. An instance myProblem is setup on a constant size time-step grid with arbitrary initial conditions and parameters.

RigidBodyProblem myProblem;

double h = 0.1;
int n_steps = 1000;
myProblem.baselinstep(0.0,h,n_steps);

Eigen::Matrix<double,3,1> Iinv(3.0/2.0,1.0,0.5);
myProblem.Iinv(Iinv);

Eigen::Matrix<double,3,1> pos(cos(M_PI/3.0), 0.0, sin(M_PI/3.0));
myProblem.pos(0,pos);

Finally an integrator is created from the RKMK::Factory template class by calling the createIntegrator function. The templates arguments should be the same than those for the RKMK::Abstract::Problem, and the function receives as arguments the instance of the problem and a string denoting the specific numerical scheme to use (in this case the classical order 4 Runge-Kutta, see REF for a list of available schemes).

RKMK::Abstract::Integrator& integrator = RKMK::Factory<double,Q,Algebra>::createIntegrator(myProblem,"RK 4");

bool success = integrator.integrate();

The result of the above code give the time evolution of the angular momentum of the body. As can clearly be seen on the following figure, the solution lies on the sphere of radius \(\Vert Y_0\Vert\), preserving the norm of the angular momentum.

Variational module