numerical integration

Numerical Integration for Rotational Dynamics

In order to advance a simulation forward in time, a numerical integration scheme is needed to evaluate the continuous part of the dynamics, which is specified by a differential equation. For translational dynamics, this is usually of the form

\begin{array}{l}\dot x = v\\\dot v=a(x,v),\end{array}

where x is the position, and v is the velocity. The acceleration, defined as \dot v, is a function of the position and velocity. Under the assumption that acceleration is constant, the differential equations have a closed-form solution, giving the update formula

\begin{array}{l}\Delta v = a \Delta t\\ \Delta x = \frac 12 (v + v') \Delta t,\end{array}

where x' = x + \Delta x and v' = v + \Delta v.

In this article, we will derive an analogous result for rotational dynamics — the numerical integration scheme for the differential equations describing rotational motion under the assumption of constant angular acceleration.

The assumption of constant acceleration is not always appropriate for translational dynamics, particularly for systems like springs that have rapidly varying acceleration. For rotational dynamics, a proper discussion of the appropriateness of constant angular acceleration is best left to a dedicated article, as it would involve the details of rigid body dynamics, which is not the concern of this article.

Translational and Rotational Analogies in 3D

Whereas translational motion describes the change in an object’s position over time, rotational motion describes the change in an object’s orientation. Positions are described relative to an origin by a displacement vector, with the null vector describing the origin. In an analogous fashion, orientations are described relative to a standard orientation by a rotation, with the identity rotation describing the standard.

Rotations in 3D can be quantified by:

  • an axis and an angle;
  • by virtue of being linear maps: a rotation matrix — a real 3-by-3 matrix with determinant 1 and inverse coinciding with its transpose;
  • a unit quaternion.

We will assume that the reader is already familiar with rotations and their various descriptions.

In 3-dimensional space, the space of rotations also possesses a 3-dimensional structure. This is not the case in general — for example, in 2D the space of rotations is 1-dimensional. Because of this, there are many analogies between translational and rotational dynamics in 3D.

Translational Rotational
position, x\in\mathbb R^3 rotation matrix, R\in\mathrm{SO}(3,\mathbb R), or unit quaternion, q\in S^3\subset\mathbb H
velocity, v=\dot x\in\mathbb R^3 angular velocity, \omega\in\mathbb R^3
acceleration, a=\dot v\in\mathbb R^3 angular acceleration, \alpha=\dot\omega\in\mathbb R^3
mass, m\in\mathbb R moment of inertia tensor, I\in\mathbb R^3\otimes\mathbb R^3\cong\mathbb R^{3\times 3}
momentum, p=mv\in\mathbb R^3 angular momentum, L=I\omega\in\mathbb R^3
force, F=\dot p\in\mathbb R^3 torque, \tau=\dot L\in\mathbb R^3

The aim of this article is not to discuss how accelerations are produced by forces, nor how angular accelerations are produced by torques, and therefore we will be focusing on the top half of this table.

The main reason that rotational motion is more complicated than translational motion is that the space of rotations is not a linear space — it is `curved’. In a linear space, a constant velocity moves you along a straight line; in the space of rotations, a constant angular velocity moves you along a curve.

The Rotational Equations

The key piece of information missing from the table above is the relation between the angular velocity, \omega, and the time derivative of the rotation, \dot R (or \dot q).

For rotation matrices this is
\dot R=\omega_{\times} R,

where \omega_{\times} is an anti-symmetric matrix describing the linear map x\mapsto\omega\times x, called the cross-product matrix,
\omega_{\times} = \left[\begin{array}{ccc}0 & -\omega_3 & \omega_2\\\omega_3 & 0 & -\omega_1\\-\omega_2 & \omega_1 & 0\end{array}\right].

For the mathematically literate, \omega_{\times} is an element of the Lie algebra \mathfrak{so}(3,\mathbb R), the tangent space to \mathrm{SO}(3,\mathbb R) at the identity. The conversion \omega\mapsto\omega_{\times} is an isomorphism between our physical space \mathbb R^3 and \mathfrak{so}(3,\mathbb R), which is possible as they are of the same dimension. Under this identification, the cross product coincides with the matrix commutator, which is the Lie bracket for the algebra. The time derivative \dot R is a tangent vector at R, which is given by the pushforward of \omega_{\times} by right translation.

The quaternionic version of this is
\dot q=\frac 12\omega q,

where \omega is interpreted as an imaginary quaternion and the multiplication is quaternionic. The identification of vectors with imaginary quaternions is analogous to the identification of vectors with anti-symmetric matrices, as the imaginary quaternions form the tangent space to the unit quaternions at the identity, 1. The factor of \frac 12 comes from the property that each rotation is described by two unique unit quaternions, whereas there is only one rotation matrix.

Integrating the Rotational Equations

The easy part of the integration is the angular velocity, which, assuming constant angular acceleration gives

\omega'=\omega + \alpha\Delta t.

However the other equation is not so straightforward.

In the special case that \alpha=0, the angular velocity is constant and the solution is given by the exponential map

R'=\exp(\omega_{\times}\Delta t) R
q'=\exp(\frac 12 \omega\Delta t) q.

But when \alpha\neq 0 and \omega varies with time, this becomes much more complicated. The solution can be written as an infinite series via the Magnus expansion,

R'=\exp( \Omega_{\times} ) R
\Omega=\Omega_1 + \Omega_2 + \cdots.

For constant angular acceleration, we can compute the first few of these terms by hand,

\Omega_1 = \frac 12 (\omega+\omega')\ \Delta t
\Omega_2 = \frac{1}{12}\ \alpha\times\omega\ \Delta t^3 = \frac{1}{12}\ \omega'\times\omega\ \Delta t^2
\Omega_3 = \frac{1}{240}\ \alpha\times(\alpha\times\omega)\ \Delta t^5.

This series is guaranteed to converge as long as the angle rotated per time-step is less than \pi/\sqrt{2} (which is about 127 degrees). Typically the first one or two terms are plenty. Note that we have given \Omega in vector form, but for use in the matrix exponential the sum should be converted to the cross-product matrix form. We can also reverse-engineer the quaternion version from this,

q'=\exp( \frac 12 \Omega ) q,

where \Omega is now an imaginary quaternion.

A Geometric Interpretation

We can think of the Magnus expansion \Omega as a polynomial in \Delta t with values in \mathfrak{so}(3,\mathbb R), the tangent space to \mathrm{SO}(3,\mathbb R) at the identity. In other words it is a path in the tangent space that begins at the null vector \Omega(0)=0 and ends at some vector \Omega(\Delta t).

The exponential map maps this tangent space to the manifold, \mathfrak{so}(3,\mathbb R)\to\mathrm{SO}(3,\mathbb R), in a way that identifies the null vector with the identity. Therefore it `wraps’ the path around the manifold so the starting point is the identity, \exp(\Omega(0))=\exp(0)=\mathrm{id}, and the end point is \exp(\Omega(\Delta t)).

The act of multiplying this path with the current rotation on the right \exp(\Omega)\mapsto \exp(\Omega)R is called right translation. It `moves’ the path across the manifold so that the starting point is now at R, and the end point is the new rotation R'. Therefore the integration moves us along paths that are constrained to rotation group.

Evaluating the Exponential Map

For both matrices and quaternions, the exponential map is defined by the series

\exp( X ) = \sum_{n=0}^\infty \frac{X^n}{n!}.

However, for the case that X is an anti-symmetric matrix, the exponential map corresponds to the conversion from axis-angle to rotation matrix (the length of the associated vector is the angle, the direction is the axis). The analogous statement applies for the quaternionic version: the exponential map acting on imaginary quaternions corresponds to the axis-angle to unit quaternion conversion,

\exp( X ) = \cos\left(\left\lVert X\right\rVert\right) + \frac{X}{\left\lVert X\right\rVert}\sin\left(\left\lVert X\right\rVert\right).

However, performing such a conversion involves evaluating trig functions, which are expensive. Instead, we can truncate the series expansion of the exponential map to get an approximation. This is appropriate as we will be evaluating it with \Delta t, so the higher powers will be negligible.

For an imaginary quaternion, X^2 = -\left\lVert X\right\rVert^2. This means the quaternion exponential map to third-order is
\exp( X ) \approx (1 - \frac 12 \left\lVert X\right\rVert^2) + (1 - \frac 16 \left\lVert X\right\rVert^2) X.

The full exponential map produces unit quaternions from imaginary quaternions, so we should also build a normalization into this approximation,
\exp( X ) \approx \frac{(1 - \frac 12 \left\lVert X\right\rVert^2) + (1 - \frac 16 \left\lVert X\right\rVert^2) X }{ \sqrt{(1 - \frac 12 \left\lVert X\right\rVert^2)^2 + (1 - \frac 16 \left\lVert X\right\rVert^2)^2 \left\lVert X\right\rVert^2}}.

This approximation is good for small rotations per time-step. The error is less than 1% for rotations up to 90 degrees per time-step. For faster rotations, the full exponential map should be applied.

Unfortunately, rotation matrices are not so easily (ortho)normalizable, and so you’re best off using the full axis-angle to rotation matrix conversion to implement the exponential map for matrices. For this reason the quaternion description is preferable to rotation matrices.

Conclusion

We have derived a method for time-stepping rotational dynamics that expresses each update as a rotation. The update rotation is determined in vector form (axis-angle) by truncating the Magnus expansion and then converted to the underlying algebraic representation (rotation matrix/unit quaternion) by applying the exponential map.

For the quaternion version, we also approximated the exponential map in a way that avoids evaluating trig functions while preserving normalization, so the resulting quaternion is always a valid rotation.

Advertisements