# Rotation matrix of Euler's equations of rotation relative to inertial reference frame

I was playing with simulation of Euler's equations of rotation in MATLAB,

$$I_1\dot{\omega}_1 + (I_3 - I_2)\omega_2\omega_3 = M_1,$$

$$I_2\dot{\omega}_2 + (I_1 - I_3)\omega_3\omega_1 = M_2,$$

$$I_3\dot{\omega}_3 + (I_2 - I_1)\omega_1\omega_2 = M_3,$$

where $I_1$, $I_2$ and $I_3$ are the principal moments of inertia of the rigid body, $\omega_1$, $\omega_2$ and $\omega_3$ are the angular velocities around the axes of these moments of inertia, $\dot{\omega}_i$ denotes the time derivative of the angular velocity $\omega_i$ and $M_i$ denotes the external torque applied along the axis of $\omega_i$.

I would also like to find the rotation of the body, but the equations above have a rotating reference frame, so finding the rotation does not seem to have a simple answer. In order to express the rotation of the body I would think that a rotation matrix as a function of time, $R(t)$, would be a good option, such that a point $\vec{x}_0$ on the rigid body can be relocated in a inertial frame of reference with,

$$\vec{x}(t) = R(t) \vec{x}.$$

This rotation matrix should change over time as the body rotates, but any two rotations can be combined into one effective rotation by multiplying the two rotation matrices. Thus the rotation matrix after a discrete time step should be,

$$R_{n+1} = D R_n,$$

where $D$ is the rotation during that time step.

Any rotation matrix (I believe with the exception of reflections) can be written as,

$$R = \begin{bmatrix} \cos\theta\!+\!n_x^2(1\!-\!\cos\theta) & n_xn_y(1\!-\!\cos\theta)\!-\!n_z\sin\theta & n_xn_z(1\!-\!\cos\theta)\!+\!n_y\sin\theta \\ n_xn_y(1\!-\!\cos\theta)\!+\!n_z\sin\theta & \cos\theta\!+\!n_y^2(1\!-\!\cos\theta) & n_yn_z(1\!-\!\cos\theta)\!-\!n_x\sin\theta \\ n_xn_z(1\!-\!\cos\theta)\!-\!n_y\sin\theta & n_yn_z(1\!-\!\cos\theta)\!+\!n_x\sin\theta & \cos\theta\!+\!n_z^2(1\!-\!\cos\theta) \end{bmatrix},$$

where $\vec{n}=\begin{bmatrix}n_x & n_y & n_z\end{bmatrix}^T$ is the axis of rotation and $\theta$ the angle of rotation.

For an infinitesimal time step the rotation, $D$, has the axis of rotation equal to the normalized angular velocity vector and the angle of rotation equal to the magnitude of the angular velocity vector times the times step. Do I need to use the angular velocity vector in the rotating or inertial reference frame for this?

I was not completely sure if I understood the notation of the answer of David Hammen, so I performed a simple test. The test involves applying two rotations. Initially the two reference frames are lined up, where x, y and z axes in the figures for the inertial reference frame always point to the left, right and up respectively, while for the rotating reference frame the x, y and z axes are represented by $\vec{e}_1$, $\vec{e}_2$ and $\vec{e}_3$ respectively. The first rotation is 90° around the x axis:

Because in both reference frames the rotation is around the x axis, thus the rotation matrices of this rotation are,

$$R_{I:0\to 1} = R_{R:0\to 1} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & -1\\ 0 & 1 & 0 \end{bmatrix},$$

where the subscripts $I$ and $R$ stand for the inertial and rotational reference frames respectively.

The next rotation is 90° around the z axis or $\vec{e}_2$:

The two reference frames now differ, thus the rotation matrices of this rotation are,

$$\begin{array}{rrrrrr} R_{I:1\to 2} = \begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix} & & & & & R_{R:1\to 2} = \begin{bmatrix} 0 & 0 & 1 \\ 0 & 1 & 0 \\ -1 & 0 & 0 \end{bmatrix} \end{array}.$$

The final orientation should look like this:

The corresponding total rotation matrix is,

$$R_{0\to 2} = \begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}.$$

When combining the two rotations in both reference frames, then the total rotation matrix can be obtained with the following order of multiplications,

$$R_{0\to 2} = R_{R:0\to 1} R_{R:1\to 2} = R_{I:1\to 2} R_{I:0\to 1},$$

thus it appears that the following is true,

$$R_{n+1} = R_n D_R = D_I R_n.$$

Since Euler's equations use the rotational reference frame, thus $D_R$ will be used. The rotation, $D_R$, using the general equation for a rotation matrix, can be simplified, by using linear approximation since $dt$ is assumed to be small, to,

$$R(t+dt) = R(t) \begin{bmatrix} 1 & -\omega_3dt & \omega_2dt \\ \omega_3dt & 1 & -\omega_1dt \\ -\omega_2dt & \omega_1dt & 1 \end{bmatrix} = R(t) \left( \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} + \begin{bmatrix} 0 & -\omega_3 & \omega_2 \\ \omega_3 & 0 & -\omega_1 \\ -\omega_2 & \omega_1 & 0 \end{bmatrix}dt \right)$$

Or in terms of the time derivative,

$$\dot{R}(t) = \frac{R(t+dt)-R(t)}{dt} = R(t) \begin{bmatrix} 0 & -\omega_3 & \omega_2 \\ \omega_3 & 0 & -\omega_1 \\ -\omega_2 & \omega_1 & 0 \end{bmatrix}.$$

Question: When numerically integrating this, together with Euler's equation of rotation, is there a way to ensure that the determinant of $R$ remains equal to one (otherwise $\vec{x}(t)$ will also be scaled)?

• If you want help understanding the Engø formula in my answer contact me on the email at my website. It look horrendous, but it is ultimately very simple and works superbly for calculating products of rotations. You shouldn't have too many problems with what you have proposed though. The Euler equns in Mathematica or MATLAB are indeed a great deal of fun. Jul 30 '15 at 7:06
• It is all explained here: cs.cmu.edu/~baraff/sigcourse/notesd1.pdf Jul 30 '15 at 13:08
• I thought you don't integrate $\dot{R}$ directly, but integrate the euler angles (and their derivatives). Jul 30 '15 at 13:15
• @ja72 - That's a recipe for massive loss of precision. While Euler angles are handy for communicating with people, there are a myriad of other representations of SO(3), all of which are better than Euler angles. Aug 1 '15 at 13:24

Do I need to use the angular velocity vector in the rotating or inertial reference frame for this?

Yes. You can do it either way.

I start with the expression that relates the time derivative of a vector quantity $\boldsymbol u$ in the inertial and rotating frames:

$$\left(\frac {d\boldsymbol u}{dt}\right)_I = \left(\frac {d\boldsymbol u}{dt}\right)_R + \boldsymbol \omega \times \boldsymbol u$$

Expressing this in coordinate spaces,

$$\dot{\boldsymbol u}_I = \mathrm T_{R\to I} \left(\dot{\boldsymbol u}_R + \boldsymbol \omega_R \times \boldsymbol u_R \right) = \mathrm T_{R\to I}\,\dot{\boldsymbol u}_R + \mathrm T_{R\to I} \mathrm{Sk}[\boldsymbol \omega_R]\,\boldsymbol u_R$$

where $\mathrm T_{R\to I}$ is the matrix that transforms a vector from its representation in the rotating frame to its corresponding representation in the inertial frame and $\mathrm{Sk}[\boldsymbol x]$ is the skew symmetric such that $\mathrm{Sk}[\boldsymbol x]\,\boldsymbol a = \boldsymbol x \times \boldsymbol a$ (in other words, the cross product matrix).

Differentiating $\boldsymbol u_I = \mathrm T_{R\to I}\,\boldsymbol u_R$ with respect to time yields

$$\dot{\boldsymbol u}_I = \mathrm T_{R\to I}\,\dot{\boldsymbol u}_R + \dot{\mathrm T} _{R\to I}\,\boldsymbol u_R$$

From which $\dot{\mathrm T} _{R\to I} = \mathrm T_{R\to I}\,\mathrm{Sk}[\boldsymbol \omega_R]$. You can similarly derive that $\dot{\mathrm T} _{I\to R} = -\mathrm T_{I \to R}\,\mathrm{Sk}[\boldsymbol \omega_I]$. Transposing yields alternate expressions. Summarizing,

\begin{aligned} \dot{\mathrm T} _{R\to I} &= \mathrm T_{R\to I}\,\mathrm{Sk}[\boldsymbol \omega_R] \\ &= \mathrm{Sk}[\boldsymbol \omega_I]\,\mathrm T_{R\to I} \\ \dot{\mathrm T} _{I \to R} &= - \mathrm T_{I \to R}\,\mathrm{Sk}[\boldsymbol \omega_I] \\ &= - \mathrm{Sk}[\boldsymbol \omega_R]\,\mathrm T_{I \to R} \end{aligned}

When numerically integrating this, together with Euler's equation of rotation, is there a way to ensure that the determinant of R remains equal to one (otherwise x⃗ (t) will also be scaled)?

Yes, using Lie group methods. Explaining this is a bit long. 140 pages long. I'll instead refer you to Iserles, Arieh, et al. "Lie-group methods." Acta Numerica 2000 9 (2000): 215-365.

• I initially had to wrap my head around how different reference frames affect the total rotation, I added a simple test case for this to my question, since I was not 100% sure if I understood your notation. Also should a symplectic integrator preserve the determinant of the rotation matrix? Aug 1 '15 at 2:05
• Might there also be a possibility to incorporate the angular momentum into this, especially when calculating the torque free solution, such that drift can be prevented? Aug 5 '15 at 12:20

To add to David Hammen's answer on the question:

When numerically integrating this, together with Euler's equation of rotation, is there a way to ensure that the determinant of $R$ remains equal to one (otherwise $\vec{x}(t)$ will also be scaled)?

Method 1 Dumb But Effective Naïve Multiplication

Whilst you are getting up to speed with more sophisticated methods, you'll get good mileage out of simply, naïve multiplying matrix exponentials of skew-symmetric matrices together in the way proposed in your question, particularly in double precision: depending on the exact precision you need, it takes a huge number of multiplications for the variance of the eigenvalues to be significant in comparison with the rotation operator in question. The variance of the numerical errors only grows linearly with iteration number (since the eigenvalues of a rotation matrix are all on the unit circle). After all, naïve implementations of the Beam Propagation Method and the Split Step Fourier Method work very well up hundreds of thousands of iterations, and they are essentially multiplying nominally unitary matrices together, so they're working in exactly the way you propose.

Method 2: Find Nearest Orthogonal Matrix

If you want to nip numerical errors in the bud, then you simply find the least squares nearest orthogonal matrix to a given matrix. For a matrix $R$ near an orthogonal matrix, this works as follows: (1) Take the matrix logarithm - MATLAB will handle this easily (or see Method 3) to map $R=e^H$ to $H$; (2) from the logarithm $H$ form $\tilde{H} = \frac{1}{2}(H - H^T)$; note that $\tilde{H}$ is skew-symmetric by construction and is the nearest skew symmetric matrix to $H$ in the least squares Frobenius norm sense (3) exponentiate back to the rotation group $SO(3)$ to find $\tilde{R} = \exp\left(\frac{1}{2}(\log R - (\log R)^T)\right)$. Now keep going with $\tilde{R}$ instead of $R$.

Method 3: Closed Form Campbell-Baker-Hausdorff Series in $SO(3)$

Here we work in the Lie algebra, noting that, for small increment rotations $\Delta R =e^Y$ applied to your overall rotation matrix $R = e^X$ where $X,\,Y$ are skew-symmetric and $Y$ is small (has small norm) there is the famous Campbell-Baker-Hausdorff series for the logarithm $Z$ of the product, i.e. $R\,\Delta R = e^X\,e^Y = e^Z$ where:

$$Z = X + Y + \frac{1}{2}[X,\,Y] +\frac{1}{12}([X,\,[X,\,Y]]+[Y,\,[Y,\,X]]) + \text{higher order bracket series of }X,\,Y$$

So we can use the CBH series to work out each step of your calculation wholly in the Lie algebra of skew-symmetric $3\times 3$ matrices.

Whilst there is no doubt that the main method presented in David Hammen's answer - i.e. the mapping of the terms of the CBH / Magnus series to rooted trees and iterating over a lexicographic ordering of trees to evaluate the most significant terms in the series - is the fastest and numerically soundest method in the general Lie group setting, in the group $SO(3)$ of rotation matrices there is a neat closed form expression for the Campbell-Baker-Hausdorff formula that is MUCH simpler. It reads:

$$Z(X,\,Y) = \frac{\arcsin(d(X,\,Y))}{d(X,\,Y)} \left(a(X,\,Y)\,\frac{X}{\|X\|}+b(X,\,Y)\,\frac{Y}{\|Y\|}+c(X,\,Y)\,\frac{X\times Y}{\|X\|\|Y\|}\right)$$

where:

$$a(X,\,Y)=\sin\|X\| \cos^2\left(\frac{\|Y\|}{2}\right)- \sin\|Y\|\sin^2\left(\frac{\|X\|}{2}\right)\frac{X\cdot Y}{\|X\|\,\|Y\|}$$ $$b(X,\,Y) = a(Y,\,X)$$ $$c(X,\,Y)=\frac{1}{2}\sin\|X\| \sin\|Y\|- 2\,\sin^2\left(\frac{\|X\|}{2}\right)\,\sin^2\left(\frac{\|Y\|}{2}\right)\frac{X\cdot Y}{\|X\|\,\|Y\|}$$ $$d(X,\,Y)=\sqrt{2\, a(X,Y)\, b(X,Y)\,\frac{X\cdot Y}{\|X\|\,\|Y\|}+a(X,Y)^2+b(X,Y)^2+c(X,Y)^2 \left(1-\frac{(X\cdot Y)^2}{\|X\|^2 \|Y\|^2}\right)}$$

and we represent the skew-symmetrix $3\times 3$ matrix $\left(\begin{array}{ccc}0&-\omega_z&\omega_y\\\omega_z&0&-\omega_x\\\omega_y&\omega_x&0\end{array}\right)$ by the column vector $\left(\begin{array}{c}\omega_x\\\omega_y\\\omega_z\end{array}\right)$ and the norm $\|\omega\|$ is the Pythagorean norm $\sqrt{\omega_x^2+\omega_y^2+\omega_z^2}$.

This formula is elegantly derived in:

K. Engø, "On the BCH-Formula in SO(3)", BIT Numerical Mathematics 41 (2001), no.3, pp629--632.

and comes from the Rodrigues formula for $SO(3)$:

$$e^H = \mathrm{id} + \frac{\sin\|H\|}{\|H\|} H + \frac{1-\cos\|H\|}{\|H\|^2}\,H^2$$

by taking heed that the skew-symmetric part of this formula is $\frac{\sin\|H\|}{\|H|} H$, whence we can find the matrix logarithm of $e^H$ by finding its skew-symmetric part $L=\frac{1}{2}(e^H -(e^H)^T)$, then calculating its magnitude (as a column vector representation) and noting that this magnitude is $\sin\|H\|$. So, the logarithm is $\arcsin\|L\| \frac{L}{\|L\|}$.

• Very nice. There appears to be a cabal of Lie algebraists in Norway and Cambridge who are developing this stuff. I used much the same as the above, but with quaternions rather than matrices. I was on the way to writing a paper when I abruptly switched fields. Trying to be a good boy, I wiped my computer clean, and it's all gone, all except for some Lorem ipsum dolor sit amet stuff sent to a coauthor early on. I did get approval to work on the paper, but its gone. Expanding the lorem ipsum to SO(3) Lie group integrators with quaternions in my spare time is not high priority. Jul 30 '15 at 9:20
• I tried to find the nearest orthogonal matrix myself and come across the polar decomposition, such that $\tilde{R} = R \left(R^T R\right)^{-\frac12}$. I did not find anywhere else a mentioning of you method with exponentials and logarithms. The results of both are similar but often slightly different. I believe the polar decomposition is better, so would there be an advantage of using the method you mentioned, maybe in computation cost? Jan 26 '16 at 16:04
• @fibonatic Actually I'm guessing they'll be very alike in computational cost: you need to find the singular value decomposition for the polar decomposition and diagonalize then find the logarithm of the eigenvalues for the nearest unitary matrix method. You'll probably need to try to answer definitively. Jan 26 '16 at 21:21
• @WetSavannaAnimalakaRodVance I did some testing in MATLAB and found that the polar decomposition is about 10 times faster than your method, but might be possible to optimize it more because I had to do a second decomposition when the exponential is used (because I did not see an obvious way that allows you to decompose the sum of the logarithms given the first decomposition). Jan 26 '16 at 23:48