Rigid Body Dynamics

Recap: Lagrangian Mechanics

Rigid body dynamics is an application of Lagrangian Mechanics to systems with specific geometric constraints.

Let the state be represented by generalized coordinates $\mathbf{q}\in Q$ and the Lagrangian be $L(\mathbf{q},\dot{\mathbf{q}})=K(\mathbf{q},\dot{\mathbf{q}})-U(\mathbf{q})$. The equations of motion are derived from the Least Action Principle, which yields the Euler-Lagrange equations:

$$\frac{\partial L}{\partial \mathbf{q}}-\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\mathbf{q}}}\right)=0$$

When applying the Lagrangian to the optimality for least action, we get

$$\underset{\text{Change of momentum}}{\underbrace{\frac{d}{dt}\left(\frac{\partial K}{\partial \dot{\mathbf{q}}}\right)}}=\underset{\text{Fictitious force}}{\underbrace{\frac{\partial K}{\partial \mathbf{q}}}}\underset{\text{Force}}{\underbrace{-\frac{\partial U}{\partial \mathbf{q}}}}$$

If there is only kinetic energy and no potential energy then the resulting motion is called pure inertial motion.

Lagrangian in Newtonian Mechanics

The generalized coordinate $\mathbf{q}$ is a low dimensional parameterization of infinitely many atoms $a\in \mathcal{A}$. There is¹ a smooth function describing the physical 3D position of each atom as a function of $\mathbf{q}$:

$$\begin{array}{rl}{\mathbb{R}}^m& \stackrel{{\mathbf{f}}_a}{\to }{\mathbb{R}}^3\\ {\mathbf{x}}_a& ={\mathbf{f}}_a(\mathbf{q})\in {\mathbb{R}}^3\end{array}$$

The kinetic and potential energy on $(\mathbf{q},\dot{\mathbf{q}})$ by pulling back the energies from 3D space

$$K(\mathbf{q},\dot{\mathbf{q}})={\int }_{a\in \mathcal{A}}\frac{1}{2}|d{\mathbf{f}}_a{|}_{\mathbf{q}}[[\dot{\mathbf{q}}]]{|}_{{\mathbb{R}}^3}^{2}daU(\mathbf{q},\dot{\mathbf{q}})={\int }_{a\in \mathcal{A}}u({\mathbf{f}}_a(\mathbf{q}))da$$

In general, the kinetic energy is a positive definite quadratic function of $\dot{\mathbf{q}}$². This defines an inner product

$$\begin{array}{rl}K(\mathbf{q},\dot{\mathbf{q}})& =\frac{1}{2}{♭}_{\mathbf{q}}(\dot{\mathbf{q}},\dot{\mathbf{q}})\\ & =\frac{1}{2}{\dot{\mathbf{q}}}^{\top }\mathbf{M}(\mathbf{q})\dot{\mathbf{q}}\end{array}$$

Definition (Momentum, Generalized Coordinate)

The momentum for the generalized coordinate $\mathbf{q}$ is defined as

$$\mathbf{p}=\frac{\partial K}{\partial \dot{\mathbf{q}}}={♭}_{\mathbf{q}}(\dot{\mathbf{q}})={\mathbf{M}}_{\mathbf{q}}\dot{\mathbf{q}}$$

This notion of momentum is the pullback of the momentum $1-$form from the moment of the classical particle system ($p=mv$).

Rigid Body Kinematics

A rigid body is a system of particles that maintains fixed distances between all pairs of particles. For a rigid body, the degrees of freedom is given by the space of Euclidean transformations:

$$Q=\mathrm{SE}(3)=\left\{\left[\begin{array}{cccc}& & & |\\ & {\mathbf{R}}_{3\times 3}& & \mathbf{c}\\ & & & |\\ 0& 0& 0& 1\end{array}\right]|{\mathbf{R}}^{\top }\mathbf{R}=I,\det (\mathbf{R})=1,\mathbf{c}\in {\mathbb{R}}^3\right\}=\{(\mathbf{R},\mathbf{c})\}$$

Here, $\mathbf{R}$ is a rotation matrix in the special orthogonal group $\mathrm{SO}(3)$ and $\mathbf{c}$ is a translation vector. $\mathrm{SE}(3)$ represents the special Euclidean group, which is the group of all rigid transformations in three-dimensional space. Each element $\mathbf{q}\in (\mathbf{R},\mathbf{c})\in Q$ is a transformation that maps a position in the body coordinate to a position in the world coordinate. For example,

$"\\usepackage{amsmath}\n\\usepackage{amssymb}\n\\usepackage{amsfonts}\n\\usepackage{xcolor}\n\n\\begin{document}\n\\begin{tikzpicture}[scale=1]\n\n$

source code

is described by the transformation

$$\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ 1\end{array}\right]=\left[\begin{array}{cccc}& & & |\\ & {\mathbf{R}}_{3\times 3}& & \mathbf{c}\\ & & & |\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\\ 1\end{array}\right]$$

To be more preicse, $Q=\mathrm{SO}(3)\times {\mathbb{R}}^3=\mathrm{SE}(3)$.

Tangent Vectors in Rigid Bodies

How can we describe a tangent vector $\dot{\mathbf{q}}\in T_{\mathbf{q}}Q$?

Differential

We can consider $(\dot{\mathbf{R}},\dot{\mathbf{c}})$ but we’d need to account for the constraint that ${\mathbf{R}}^{\top }\mathbf{R}=I$.

Recall that since $\mathbf{R}$ is a rotation matrix, it must be orthogonal, meaning that ${\mathbf{R}}^{\top }\mathbf{R}=I$. So, taking the derivative of both sides with respect to time, we get

$$\begin{array}{rl}\frac{d}{dt}({\mathbf{R}}^{\top }\mathbf{R})& =\frac{d}{dt}I\\ {\dot{\mathbf{R}}}^{\top }\mathbf{R}+{\mathbf{R}}^{\top }\dot{\mathbf{R}}& =0\\ ({\mathbf{R}}^{\top }\dot{\mathbf{R}}{)}^{\top }+{\mathbf{R}}^{\top }\dot{\mathbf{R}}& =0\end{array}$$

Therefore ${\mathbf{R}}^{\top }\dot{\mathbf{R}}$ is a skew-symmetric matrix. Similarly,

$${\mathbf{R}}^{\top }\mathbf{R}=I⟹\mathbf{R}{\mathbf{R}}^{\top }=I$$

Taking the derivative of both sides with respect to time, we get

$$\begin{array}{rl}\frac{d}{dt}(\mathbf{R}{\mathbf{R}}^{\top })& =\frac{d}{dt}I\\ \dot{\mathbf{R}}{\mathbf{R}}^{\top }+\mathbf{R}{\dot{\mathbf{R}}}^{\top }& =0\\ \dot{\mathbf{R}}{\mathbf{R}}^{\top }+(\dot{\mathbf{R}}{\mathbf{R}}^{\top }{)}^{\top }& =0\end{array}$$

and so both ${\mathbf{R}}^{\top }\dot{\mathbf{R}}$ and $\dot{\mathbf{R}}{\mathbf{R}}^{\top }$ are skew-symmetric matrices. Let

$$\begin{array}{rl}\mathbf{A}& ={\mathbf{R}}^{\top }\dot{\mathbf{R}}\\ \mathbf{W}& =\dot{\mathbf{R}}{\mathbf{R}}^{\top }\end{array}$$

Then $\dot{\mathbf{R}}=\mathbf{RA}=\mathbf{WR}$. By theorem, every skew-symmetric matrix is a cross product with a vector, such that

$$\dot{\mathbf{R}}=\mathbf{R}[\Omega \times ]=[\omega \times ]\mathbf{R}$$

where

$$\mathbf{A}=[\Omega \times ]=\left[\begin{array}{ccc}0& -{\Omega }_3& {\Omega }_2\\ {\Omega }_3& 0& -{\Omega }_1\\ -{\Omega }_2& {\Omega }_1& 0\end{array}\right]\mathbf{W}=[\omega \times ]=\left[\begin{array}{ccc}0& -{\omega }_3& {\omega }_2\\ {\omega }_3& 0& -{\omega }_1\\ -{\omega }_2& {\omega }_1& 0\end{array}\right]$$

We call $\Omega \in {\mathbb{R}}_{\text{body}}^3$ the body coordinate angular velocity and $\omega \in {\mathbb{R}}_{\text{world}}^3$ the world coordinate angular velocity. Indeed, one can verify that $\omega =\mathbf{R}\Omega$.

Thus, a tangent vector $\dot{\mathbf{q}}\in T_{\mathbf{q}}Q$ is represented by a 6D vector

$$(\omega ,\dot{\mathbf{c}})\in {\mathbb{R}}_{\text{world}}^3\times {\mathbb{R}}_{\text{world}}^3(\Omega ,\dot{\mathbf{c}})\in {\mathbb{R}}_{\text{body}}^3\times {\mathbb{R}}_{\text{world}}^3$$

We can pick if we want to represent the angular velocity in the world coordinate or the body coordinate. In particular, suppose the $i-$th point in the body has position and mass

$${\mathbf{a}}_i\in {\mathbb{R}}_{\text{body}}^3{m}_i\in {\mathbb{R}}_{\ge 0}$$

then the translated position of the $i-$th point in the world coordinate ${\mathbb{R}}_{\text{world}}^3$ is

$${\mathbf{x}}_i=\mathbf{R}{\mathbf{a}}_i+\mathbf{c}$$

and the velocity of the $i-$th point in the world coordinate is

$$\begin{array}{rl}{\dot{\mathbf{x}}}_i& =[\omega \times ]{\mathbf{x}}_i+\dot{\mathbf{c}}\\ & =\mathbf{R}[\Omega \times ]{\mathbf{a}}_i+\dot{\mathbf{c}}\\ & =\mathbf{R}(\Omega \times {\mathbf{a}}_i)+\dot{\mathbf{c}}\end{array}$$

This makes sense because we get the world angular velocity contribution from $[\Omega \times ]$. Then we can convert to the body angular velocity contribution by multiplying with $\mathbf{R}$, and finally we apply theorem to show the cross product.

With out velocity, we can compute the kinetic energy as

$$\begin{array}{rl}K(\mathbf{x},\dot{\mathbf{x}})& =\sum _{i\in \text{body}}\frac{m_i}{2}|{\dot{\mathbf{x}}}_i{|}^2\\ & =\sum _{i\in \text{body}}\frac{m_i}{2}{\left|\mathbf{R}(\Omega \times {\mathbf{a}}_i)+\dot{\mathbf{c}}\right|}^2\\ & =\sum _{i\in \text{body}}\frac{m_i}{2}{\left|\mathbf{R}(\Omega \times {\mathbf{a}}_i)\right|}^2+\sum _{i\in \text{body}}{m}_i⟨\mathbf{R}(\Omega \times {\mathbf{a}}_i),\dot{\mathbf{c}}⟩+\sum _{i\in \text{body}}\frac{m_i}{2}|\dot{\mathbf{c}}{|}^2\\ & =\sum _{i\in \text{body}}\frac{m_i}{2}{\left|\Omega \times {\mathbf{a}}_i\right|}^2+\sum _{i\in \text{body}}{m}_i⟨\mathbf{R}(\Omega \times {\mathbf{a}}_i),\dot{\mathbf{c}}⟩+\frac{m}{2}|\dot{\mathbf{c}}{|}^2\end{array}$$

In the first term, we can remote $\mathbf{R}$ because the norm of a rotation matrix is $1$ (it preserves lengths). In the third term, $m$ is simply the total mass of the body.

Consider the second term. All the operations listed are linear, and that $\mathbf{R},\Omega$ are not dependent on particle $i$, so we can factor them out of the summation:

$$\sum _{i\in \text{body}}{m}_i⟨\mathbf{R}(\Omega \times {\mathbf{a}}_i),\dot{\mathbf{c}}⟩=⟨\mathbf{R}\left(\Omega \times \sum _{i\in \text{body}}{m}_i{\mathbf{a}}_i\right),\dot{\mathbf{c}}⟩$$

Center of Mass

Recall that the formula for calculating the position vector of the center of mass is

$${\mathbf{r}}_{\text{com}}=\frac{1}{m}\sum _{i\in \text{body}}{m}_i{\mathbf{a}}_i$$

We can get

$$m{\mathbf{r}}_{\text{com}}=\sum _{i\in \text{body}}{m}_i{\mathbf{a}}_i$$

If we choose the body coordinate $\mathbf{a}$ such that the center of mass is at the origin, then ${\mathbf{r}}_{\text{com}}=0$ (recall it is relative to the body coordinate) and so

$$\sum _{i\in \text{body}}{m}_i{\mathbf{a}}_i=m{\mathbf{r}}_{\text{com}}=0$$

Returning back to the kinetic energy, we get

$$K(\mathbf{x},\dot{\mathbf{x}})=\sum _{i\in \text{body}}\frac{m_i}{2}{\left|\Omega \times {\mathbf{a}}_i\right|}^2+\frac{m}{2}|\dot{\mathbf{c}}{|}^2$$

Theorem (Kinetic Energy Decoupling)

If the origin of the body coordinate is at the center of mass, then the kinetic energy is

$$K(\mathbf{q},\dot{\mathbf{q}})=\underset{K_{\text{com}}(\dot{\mathbf{c}})}{\underbrace{\frac{m}{2}|\dot{\mathbf{c}}{|}^2}}+\underset{K_{\text{rot}}(\mathbf{R},\dot{\mathbf{R}})}{\underbrace{\sum _{i\in \text{body}}\frac{m_i}{2}{\left|\Omega \times {\mathbf{a}}_i\right|}^2}}$$

That is, the translation and rotation contributions to the kinetic energy are decoupled. Indeed,

• the pure translational kinetic energy of the center of mass moving through space is the first term
• the pure rotational kinetic energy about the center of mass is the second term

Definition (Inertia Tensor)

We can further expand the $K_{\text{rot}}(\mathbf{R},\dot{\mathbf{R}})$ term.

$$\begin{array}{rl}{K}_{\text{rot}}(\mathbf{R},\dot{\mathbf{R}})& =\sum _{i\in \text{body}}\frac{m_i}{2}{\left|\Omega \times {\mathbf{a}}_i\right|}^2\\ & =\sum _{i\in \text{body}}\frac{m_i}{2}{\left|[{\mathbf{a}}_i\times ]\Omega \right|}^2\\ & =\sum _{i\in \text{body}}\frac{m_i}{2}\left({\Omega }^{\top }[{\mathbf{a}}_i\times {]}^{\top }[{\mathbf{a}}_i\times ]\Omega \right)\\ & =\frac{1}{2}{\Omega }^{\top }\left(\sum _{i\in \text{body}}{m}_i[{\mathbf{a}}_i\times {]}^{\top }[{\mathbf{a}}_i\times ]\right)\Omega \\ & =\frac{1}{2}{\Omega }^{\top }\left(\sum _{i\in \text{body}}{m}_i\left(|{\mathbf{a}}_i{|}^{2}I-{\mathbf{a}}_i{\mathbf{a}}_i^{\top }\right)\right)\end{array}$$

There are many steps.

1. By the anti-commutative property of the cross product and theorem
2. Since $[{\mathbf{a}}_i\times ]\Omega$ is a vector, its squared norm is equivalent to dot product with itself.
3. We can factor out body angular velocity $\Omega$ since it is identical to every particle.
4. By theorem, we get the final line.

For a rigid body centered at the origin of the body coordinate, the moment of intertia or the inertia tensor is defined as

$${\mathbf{I}}_{\text{body}}=\sum _{i\in \text{body}}{m}_i\left(|{\mathbf{a}}_i{|}^{2}I-{\mathbf{a}}_i{\mathbf{a}}_i^{\top }\right)$$

For a continuous object, we can replace the summation with an integral:

$${\mathbf{I}}_{\text{body}}={\iiint }_{\mathbf{a}\in \text{body}}\left(|\mathbf{a}{|}^{2}I-\mathbf{a}{\mathbf{a}}^{\top }\right)d{m}_{\mathbf{a}}$$

Inertia Tensor in Kinetic Energy

With the inertia tensor, we can express the total inertial kinetic energy (assuming the center of mass is at the origin) from Theorem (Kinetic Energy Decoupling) as

$$K(\mathbf{q},\dot{\mathbf{q}})=\frac{m}{2}|\dot{\mathbf{c}}{|}^2+\frac{1}{2}{\Omega }^{\top }{\mathbf{I}}_{\text{body}}\Omega$$

where the $3\times 3$ symmetric matrix ${\mathbf{I}}_{\text{body}}$ is the moment of inertia which is time independent.

Using the world coordinate, we can express the total inertial kinetic energy as

$$K(\mathbf{q},\dot{\mathbf{q}})=\frac{m}{2}|\dot{\mathbf{c}}{|}^2+\frac{1}{2}{\omega }^{\top }{\mathbf{I}}_{\text{world}}\omega$$

where

$${\mathbf{I}}_{\text{world}}=\mathbf{R}{\mathbf{I}}_{\text{body}}{\mathbf{R}}^{\top }$$

is the inertia in the world coordinate. For a rigid body, ${\mathbf{I}}_{\text{body}}$ is time-independent. However, ${\mathbf{I}}_{\text{world}}$ is dynamic, and co-rotates with the rigid body.

Inertia Tensor in World Coordinate

In particular, we say it is dynamic because its axes of rotation changes with the rigid body’s rotation with respect to the world coordinate axes. The fact that

$${\mathbf{I}}_{\text{world}}=\mathbf{R}{\mathbf{I}}_{\text{body}}{\mathbf{R}}^{\top }$$

means that finding its position in the world coordinate, we must update ${\mathbf{I}}_{\text{world}}$ to reflect the object’s current orientation.

As the rigid body rotates, $\mathbf{R}$ changes, causing ${\mathbf{I}}_{\text{world}}$ to change accordingly. This reflects how the distribution of mass relative to the world axes changes as the body rotates, even though the distribution of mass relative to the body axes (captured by ${\mathbf{I}}_{\text{body}}$) remains constant.

Theorem (Moment of Inertia Semi-Definiteness)

The moment of inertia is positive semi-definite, and its $3$ eigenvalues satisfy the triangle inequality.

Proof:

WLOG, we can choose³ the basis that diagonalizes the inertia tensor, such that

$$\mathbf{I}=\left[\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_2& 0\\ 0& 0& {\lambda }_3\end{array}\right]$$ $$i_x=\int \frac{|x{|}^2}{2}dm{i}_y=\int \frac{|y{|}^2}{2}dm{i}_z=\int \frac{|z{|}^2}{2}dm$$

These variables will act as “intermediates” to help us express the eigenvalues of the inertia tensor.

The diagonal entries (and thus the eigenvalues) of the inertia tensor are

$${\lambda }_1={\int }_{\text{body}}\left(|y{|}^2+|z{|}^2\right)dm=i_y+i_z$$

and similarly ${\lambda }_2=i_x+i_z$ and ${\lambda }_3=i_x+i_y$. Because $i_x,i_y,i_z$ are all non-negative, then any sum of them must also be positive. Therefore ${\lambda }_1,{\lambda }_2,{\lambda }_3$ are non-negative. This proves that the inertia tensor is positive semi-definite.

Moreover, we can check that the eigenvalues satisfy the triangle inequality.

$$\begin{array}{rl}{\lambda }_1+{\lambda }_2& \ge {\lambda }_3\\ (i_y+i_z)+(i_x+i_z)& \ge i_x+i_y\\ 2i_z& \ge 0\end{array}$$

and so we are done.

Main Idea

Recall that the moment of inertia is defined as

$$I=m{r}^2$$

where $m$ is the mass and $r$ is the distance from the axis of rotation (radius). The idea is that on each axis, the moment of inertia is the sum of the contributions from all mass elements, and each contribution is proportional to the square of the distance from that axis.

In particular, if we consider the $x$-axis, then the distance from the $x$-axis for a mass element at position $(x,y,z)$ is $r^2=y^2+z^2$.

Definition (Spinning Top)

A rigid body is called a spinning top if a pair of eigenvalues of its inertia tensor are equal. That is, for the spinning top, there is a continuous family of rotation $\mathbf{Q}$ so that

$$\mathbf{Q}{\mathbf{I}}_{\text{body}}{\mathbf{Q}}^{\top }={\mathbf{I}}_{\text{body}}$$

Equations of Motion for Rigid Bodies

The equations of motion for a rigid body can be described as follows. Let

• $\mathbf{R}$ be the rotation matrix representing the orientation of the rigid body
• $\mathbf{c}$ be the translation vector representing the position of the center of mass of the rigid body
• $\omega$ be the world coordinate angular velocity
• $\Omega$ be the body coordinate angular velocity
• ${\mathbf{I}}_{\text{body}}$ be the inertia tensor in the body coordinate
• ${\mathbf{I}}_{\text{world}}$ be the inertia tensor in the world coordinate

In total,

$$\text{Body-Frame Formulation}:\{\begin{array}{ll}\dot{\mathbf{R}}& =\mathbf{R}[\Omega \times ]\\ {\mathbf{I}}_{\text{body}}\dot{\Omega }& =-\Omega \times ({\mathbf{I}}_{\text{body}}\Omega )\\ m\ddot{\mathbf{c}}& =-\frac{\partial U}{\partial \mathbf{c}}\end{array}$$

and likewise,

$$\text{World-Frame Formulation}:\{\begin{array}{ll}\dot{\mathbf{R}}& =[\omega \times ]\mathbf{R}\\ {\mathbf{I}}_{\text{world}}\dot{\omega }& =-\omega \times ({\mathbf{I}}_{\text{world}}\omega )\\ {\mathbf{I}}_{\text{world}}& =\mathbf{R}{\mathbf{I}}_{\text{body}}{\mathbf{R}}^{\top }\\ m\ddot{\mathbf{c}}& =-\frac{\partial U}{\partial \mathbf{c}}\end{array}$$

Applying the Lagrangian to Rigid Bodies

Recall that the Lagrangian is defined as $L=K-U$. If we assume there are no external forces, then the potential energy $U$ is $0$ and so $L=K$. From Inertia Tensor in Kinetic Energy, we can take the variation from only the rotational kinetic energy.

$${\int }_0^T\frac{1}{2}\Omega (t{)}^{\top }{\mathbf{I}}_{\text{body}}\Omega (t)dt$$

Theorem (Angular Momentum is Conserved in World Coordinate)

The conservation of angular momentum depends on which coordinate system we use. The world coordinate angular momentum is defined as

$$\mathbf{l}:={\mathbf{I}}_{\text{world}}\omega$$

We say that the world coordinate angular momentum is conserved if $\mathbf{l}$ is constant over time. In particular,

$$\frac{d}{dt}\mathbf{l}=0$$

Since both ${\mathbf{I}}_{\text{world}}$ and $\omega$ are dependent on time, we can apply the product rule to get

$$\begin{array}{rl}0& =\frac{d}{dt}{\mathbf{I}}_{\text{world}}\omega \\ & ={\dot{\mathbf{I}}}_{\text{world}}\omega +{\mathbf{I}}_{\text{world}}\dot{\omega }\end{array}$$

The goal is to manipulate the first term so that is equal to one of the World-Frame Formulation equations of motion. Indeed, we can compute the product rule again

$$\begin{array}{rl}{\dot{\mathbf{I}}}_{\text{world}}& =\frac{d}{dt}(\mathbf{R}{\mathbf{I}}_{\text{body}}{\mathbf{R}}^{\top })\\ & =\dot{\mathbf{R}}{\mathbf{I}}_{\text{body}}{\mathbf{R}}^{\top }+\mathbf{R}{\mathbf{I}}_{\text{body}}{\dot{\mathbf{R}}}^{\top }\end{array}$$

Recall that ${\mathbf{I}}_{\text{body}}$ is time independent, so its derivative is $0$. By equation 1 of the World-Frame Formulation, we can substitute $\dot{\mathbf{R}}=[\omega \times ]\mathbf{R}$ to get

$$\begin{array}{rl}\dot{\mathbf{R}}{\mathbf{I}}_{\text{body}}{\mathbf{R}}^{\top }+\mathbf{R}{\mathbf{I}}_{\text{body}}{\dot{\mathbf{R}}}^{\top }& =([\omega \times ]\mathbf{R}){\mathbf{I}}_{\text{body}}{\mathbf{R}}^{\top }+\mathbf{R}{\mathbf{I}}_{\text{body}}({\mathbf{R}}^{\top }[-\omega \times ])\\ & =[\omega \times ]{\mathbf{I}}_{\text{world}}-{\mathbf{I}}_{\text{world}}[\omega \times ]\end{array}$$

Thus

$${\dot{\mathbf{I}}}_{\text{world}}\omega =[\omega \times ]{\mathbf{I}}_{\text{world}}\omega -{\mathbf{I}}_{\text{world}}[\omega \times ]\omega$$

Since $[\omega \times ]\omega$ is equivalent to $\omega \times \omega$ (via theorem), and $\omega \times \omega =0$ (by anti-commutative property of the cross product), we get

$${\dot{\mathbf{I}}}_{\text{world}}\omega =[\omega \times ]{\mathbf{I}}_{\text{world}}\omega$$

Substituting back into the product rule, we get

$$0=[\omega \times ]{\mathbf{I}}_{\text{world}}\omega +{\mathbf{I}}_{\text{world}}\dot{\omega }$$

which is exactly the second equation of the World-Frame Formulation. Therefore using the world coordinate angular momentum conserves angular momentum.

Theorem (Angular Momentum is Not Conserved in Body Coordinate)

Imagine we are an ant glued to the rigid body flying around in empty space. Since we have the world frame angular momentum, translating it to the body frame is just a rotation. The angular momentum in the body frame is defined as

$$\mathbf{L}:={\mathbf{I}}_{\text{body}}\Omega$$

Recall this is similar to the classical definition of angular momentum, $L=I\omega$ where $I$ is the moment of inertia and $\omega$ is the angular velocity. From how we derived the integrand, we know kinetic energy wrt to the body is

$$K=\frac{1}{2}{\Omega }^{\top }{\mathbf{I}}_{\text{body}}\Omega$$

Thus,

$$\begin{array}{rl}\frac{dK}{d\Omega }& ={\mathbf{I}}_{\text{body}}\Omega \\ & =\mathbf{L}\end{array}$$

We can reuse from before that

$$\begin{array}{rl}\mathbf{l}& ={\mathbf{I}}_{\text{world}}\omega \\ & =(\mathbf{R}{\mathbf{I}}_{\text{body}}{\mathbf{R}}^{\top })(\mathbf{R}\Omega )\\ & =\mathbf{R}{\mathbf{I}}_{\text{body}}\Omega \\ & =\mathbf{RL}\end{array}$$

and so $\mathbf{L}={\mathbf{R}}^{\top }\mathbf{l}$. And since $\mathbf{R}$ (and thus ${\mathbf{R}}^{\top }$) is dependent on time, angular momentum is not conserved.

Numerical Methods for Rigid Body Simulation

We’ll use the World-Frame Formulation equations of motion to simulate the rigid body dynamics.

Current methods such as forward Euler or RK4 do not respect rotation matrices⁴.

Simple First-Order Method

Here is an example of a simple first-order method based on conservation of angular momentum.

Input: Initial $\mathbf{R},\omega ,{\mathbf{I}}_{\text{body}},\Delta t$

1. Compute world-space moment of inertia: ${\mathbf{I}}_{\text{world}}=\mathbf{R}{\mathbf{I}}_{\text{body}}{\mathbf{R}}^{\top }$
2. Compute world-space angular momentum: $\mathbf{L}={\mathbf{I}}_{\text{world}}\omega$ (this is conserved if no external torques)
3. For each frame: a. Update angular velocity: $\omega \leftarrow {\mathbf{I}}_{\text{world}}^{-1}\mathbf{L}$ b. Update rotation: $\mathbf{R}\leftarrow \text{Rotation}(\Delta t|\omega |,\frac{\omega }{|\omega |})\mathbf{R}$ (angle-axis rotation) c. Recompute world-space moment of inertia: ${\mathbf{I}}_{\text{world}}=\mathbf{R}{\mathbf{I}}_{\text{body}}{\mathbf{R}}^{\top }$ d. Render the object using $\mathbf{R}$ as the orientation matrix.
4. End

Buss’ Augmented Second-Order Method

From Samuel Buss 2001 "Accurate and Efficient Simulation of Rigid Body Dynamics", (link to paper). This method provides better accuracy for large time steps.

Input: Initial $\mathbf{R},\omega ,{\mathbf{I}}_{\text{body}},\Delta t$

1. Compute world-space moment of inertia: ${\mathbf{I}}_{\text{world}}=\mathbf{R}{\mathbf{I}}_{\text{body}}{\mathbf{R}}^{\top }$
2. Compute world-space angular momentum: $\mathbf{L}={\mathbf{I}}_{\text{world}}\omega$ (conserved)
3. For each frame: a. Update angular velocity: $\omega \leftarrow {\mathbf{I}}_{\text{world}}^{-1}\mathbf{L}$ b. Compute angular acceleration term: $\alpha =-{\mathbf{I}}_{\text{world}}^{-1}(\omega \times \mathbf{L})$ c. Compute effective angular velocity: $\bar{\omega }=\omega +\frac{\Delta t}{2}\alpha +\frac{(\Delta t{)}^2}{12}(\alpha \times \omega )$ d. Update rotation: $\mathbf{R}\leftarrow \text{Rotation}(\Delta t|\bar{\omega }|,\frac{\bar{\omega }}{|\bar{\omega }|})\mathbf{R}$ e. Recompute world-space moment of inertia: ${\mathbf{I}}_{\text{world}}=\mathbf{R}{\mathbf{I}}_{\text{body}}{\mathbf{R}}^{\top }$ f. Render the object using $\mathbf{R}$ as the orientation matrix.
4. End

The idea is that we use Dynkin’s formula and Baker-Campbell-Hausdorff formula to derive a more accurate effective angular velocity $\bar{\omega }$ that accounts for the non-commutativity of rotations (i.e. $e^{\mathbf{A}}{e}^{\mathbf{B}}\ne e^{\mathbf{A}+\mathbf{B}}$). Instead,

$$e^{\mathbf{A}}{e}^{\mathbf{B}}=\exp \left(\mathbf{A}+\mathbf{B}+\frac{1}{2}[\mathbf{A},\mathbf{B}]+\frac{1}{12}[\mathbf{A},[\mathbf{A},\mathbf{B}]]+\frac{1}{12}[\mathbf{B},[\mathbf{B},\mathbf{A}]]+\cdots \right)$$

where $[\mathbf{A},\mathbf{B}]=\mathbf{AB}-\mathbf{BA}$. Indeed, this means that

$$\left[[a\times ],[b\times ]\right]=\left[(a\times b)\times \right]$$

giving us line (3c) of the algorithm.

Idea

This is similar to RK2.

Footnotes

1. Why? ↩

2. Kinetic energy cannot be negative. Since kinetic energy here is expressed as a quadratic form in generalized velocities $\dot{\mathbf{q}}$, the matrix $\mathbf{M}(\mathbf{q})$ must be positive definite to ensure that $K$ is always positive for any nonzero $\dot{\mathbf{q}}$. Especially if we consider when $\dot{\mathbf{q}}=0$, then $K(\mathbf{q},0)=0$ for all $\mathbf{q}$. ↩

3. TODO: show that we can always choose such a basis. Let ↩

4. Why? ↩