$$\newcommand{\u}[1]{\boldsymbol{\mathsf{#1}}} \renewcommand{\b}[1]{\boldsymbol{#1}} \newcommand{\t}[1]{\textsf{#1}} \newcommand{\m}[1]{\mathbb{#1}} \def\RR{\bf R} \def\bold#1{\bf #1} \def\mbf#1{\mathbf #1} \def\uv#1{\hat{\usf {#1}}} \def\dl#1{\underline{\underline{#1}}} \newcommand{\usf}[1]{\boldsymbol{\mathsf{#1}}} \def\bs#1{\usf #1}$$

Dynamic Euler Equation: Lie-Group Integration Scheme

The dydnamic Euler equation is $\begin{align} \bs{J}\dot{\bs{\omega}} + \bs{\Omega}\bs{J}\omega - \bs{T} = 0, \end{align}$ where $\bs{\omega}$ and $\bs{J}$ are the rigid body angular velocity vector and inertia tensor with respect to the center of mass, expressed in the reference frame attached to the body; $\bs{T}$ is the external torque related to the body center of mass, expressed in the same reference frame; and $\begin{align} \bs{\Omega} = \dot{\bs{R}}\bs{R}^{\text{T}} \end{align}$ is the angular velocity tensor, which is a skew symmetric tensor with an axial vector, i.e., the angular velocity vector $\bs{\omega}$.

We use the integration scheme presented by Terze et al. [1] to compute the instant rotation tensor $\bs{R}$ and angular velocity vector $\bs{\omega}$ at time $t$ with the initial conditions $\bs{R}(0)=\bs{R}_0$ and $\bs{\omega}(0)=\bs{\omega}_0$.

Suppose we know the rotation tensor $\bs{R}$ and angular velocity vector $\bs{\omega}$ at $n_{th}$ time step, denoted as $\bs{R}_n$ and $\bs{\omega}_n$. Then we can update the rotation tensor and velocity vector as follows:

$\begin{align} \bs{Y}_{n+\frac{1}{2}} &= \exp\left( -\frac{h}{2} \bs{\Omega}_n \right) \left( \bs{J} \bs{\omega}_n + \frac{h}{2}\bs{T}_n \right)\\ \bs{\omega}_{n+\frac{1}{2}} &= \bs{J}^{-1} \bs{Y}_{n+\frac{1}{2}} \\ \bs{R}_{n+1} &= \bs{R}_n \exp\left(h \bs{\Omega}_{n+\frac{1}{2}}\right) \\ \bs{Y}_{n+1} &= \exp\left( -\frac{h}{2} \bs{\Omega}_{n+\frac{1}{2}} \right) \left[\exp\left( -\frac{h}{2} \bs{\Omega}_{n+\frac{1}{2}} \right) \left( \bs{J} \bs{\omega}_n + \frac{h}{2}\bs{T}_n \right) + \frac{h}{2} \exp\left( \frac{h}{2} \bs{\Omega}_{n+\frac{1}{2}} \right) \bs{T}_{n+1} \right]\\ \bs{\omega}_{n+1} &= \bs{J}^{-1} \bs{Y}_{n+1} \end{align}$ If the body in motion is free of external torque, then $\bs{T} = 0$. The $\bs{\Omega}_n$ is the skew symmetric tensor with its axial vector given by $\bs{\omega}_n$. The $\exp(\cdot)$ is the exponential map of a skew-symmetric tensor [2]. For example, given a skew symmetric tensor $\bs{B}$ with matrix representation $\begin{align} [\bs{B}] = \begin{bmatrix} 0 & -c & b \\ c & 0 & -a \\ -b & a & 0 \end{bmatrix}, \end{align}$ we have $\begin{align} \exp(\bs{B}) = \bs{I} + \frac{\sin\alpha}{\alpha} \bs{B} + \frac{1-\cos\alpha}{\alpha^2}\bs{B}^2, \end{align}$ where $\alpha = \sqrt{a^2+b^2+c^2}$.

Numerical Example

• A Tossed Book in the air

  Mathematica Code is available here

  

  Video

• A Stumbling T-handle in Space

  

  Video

• Anti-Gravity Wheel (Gyroscopic Precession)

  

  Video

References:

[1]. Zdravko Terze, Andreas Müller and Dario Zlatar, An Angular Momentum and Energy Conserving Lie-Group Integration Scheme for Rigid Body Rotational Dynamics Originating From Stormer–Verlet Algorithm, J. Comput. Nonlinear Dynam 10(5):051005, 2015.

[2]. Jean Gallier and Dianna Xu, Computing Exponentials of Skew Symmetric Matrices And Logarithms of Orthogonal Matrices, International Journal of Robotics and Automation, 17(4):1-11, 2002.

[edit]