Is kinematics required to solve for the quaternion dynamics of a tumbling body?

Question

I'm trying to solve for the rotational motion of a rigid body in the absence of external torques using quaternions in MATLAB. Assuming the axis of rotation as the unit vector $\begin{bmatrix}a_x & a_y & a_z\end{bmatrix}$ and angle of rotation as $\theta$, the quaternion variable is given by

$$\mathbf{q} = \begin{bmatrix}\cos\frac{\theta}{2} & a_x\sin\frac{\theta}{2} & a_y\sin\frac{\theta}{2} & a_z\sin\frac{\theta}{2}\end{bmatrix}^T = \begin{bmatrix}q_0 & q_1 & q_2 & q_3\end{bmatrix}^T$$

The single and double derivatives are denoted by $\mathbf{\dot{q}}$ and $\mathbf{\ddot{q}}$, respectively. In order to obtain $\mathbf{q}$ from $\mathbf{\ddot{q}}$, can I simply do a double integral to obtain the inertial states or do I need the kinematics? If the kinematics is required, then how should I incorporate it into the solver?

Dynamics equation: \begin{equation} \mathbf{\ddot{q}}=-2(\mathbf{H^TJH})^{-1}\mathbf{\dot{H}JH\dot{q}} \tag{1} \end{equation} where $\mathbf{H} = \begin{bmatrix}-2q_1 & 2q_0 & 2q_3 & -2q_2 \\ -2q_2 & -2q_3 & 2q_0 & 2q_1 \\ -2q_3 & 2q_2 & -2q_1 & 2q_0 \end{bmatrix}$ relates $\mathbf{\dot{q}}$ to the angular velocity by $\boldsymbol{\omega}=\mathbf{H\dot{q}} = \begin{bmatrix}\omega_x & \omega_y & \omega_z \end{bmatrix}^T$ and $\mathbf{J}$ is the inertia matrix.

Kinematics equation: \begin{equation} \mathbf{\dot{q}}=\displaystyle\frac{1}{2}\mathbf{\Omega q} \tag{2} \end{equation} where $\mathbf{\Omega} = \begin{bmatrix}0 & -\omega_x & -\omega_y & -\omega_z\\ \omega_x & 0 & \omega_z & -\omega_y\\ \omega_y & -\omega_z & 0 & -\omega_x \\ \omega_z & \omega_y & -\omega_x & 0\end{bmatrix}$.

Answer 1

Yes, you need the kinematical equation. What you have is a coupled system of ordinary differential equations (ODEs). ODE integrators like you'll find in MATLAB (or Mathematica, SciPy, Maple, GSL, etc.) can integrate systems of equations that look like this: \begin{align} \dot{s}_1 &= f_1(t, s_1, s_2, \ldots, s_n), \\ \dot{s}_2 &= f_2(t, s_1, s_2, \ldots, s_n), \\ \vdots \\ \dot{s}_n &= f_n(t, s_1, s_2, \ldots, s_n), \end{align} where the $s_i$ variables represent the state of your system, and you have to provide the functions $f_i$. Note that there is precisely one derivative on the left-hand side, and precisely zero derivatives on the right-hand side — although the $s_i$ may themselves represent first derivatives. But the crucial point is that your ODE integrator will only return the functions $s_1(t), \ldots, s_n(t)$; it will not return the integrals of those functions. Now, in your case, \begin{align} s_1 &= \dot{q}_0, \\ s_2 &= \dot{q}_1, \\ s_3 &= \dot{q}_2, \\ s_4 &= \dot{q}_3, \end{align} and your dynamical equation gives you expressions for $f_1, \ldots, f_4$. But plugging this system into an ODE integrator would only give you results for these quantites: $s_1(t), \ldots, s_4(t)$, which is just $\dot{\mathbf{q}}(t)$. You will not get $\mathbf{q}(t)$ itself. Now, at this point, you might reasonably think to yourself "Well, that's fine; I'll just take $\dot{\mathbf{q}}(t)$, then go back and integrate it as usual to get $\mathbf{q}(t)$." But you needed to know $\mathbf{q}(t)$ to integrate the equations in the first place, since it appears via $\mathbf{H}$ on the right-hand side of the equation for $\ddot{\mathbf{q}}$. That is, you need to evolve $\mathbf{q}(t)$ right along with the rest of the system, which means that you have to include it in your list of $s_i$'s: \begin{align} s_5 &= {q}_0, \\ s_6 &= {q}_1, \\ s_7 &= {q}_2, \\ s_8 &= {q}_3. \end{align} And the functions $f_5, \ldots, f_8$ are only given by the kinematical equation; not the dynamical.

There's an important subtlety here: Your $f_i$ equations really have to be written in terms of the $s_i$ variables; you can't have them in terms of $\dot{s}_i$ or $\dot{q}_i$. Only $s_1, \ldots, s_8$ are allowed on the right-hand sides. This may seem a bit confusing and weird, since you and I know that $s_1$ and $\dot{s}_5$, for example, are supposed to be the same things. But the ODE solver can't understand this fact, and has to treat them as completely separate and independent things.

Ultimately the reason is that we have theorems establishing the existence and uniqueness of solutions to systems of ODEs in this form (if we also have initial values). That is, the theorems tell us that the ODE integrator should work. There are no such theorems for systems like $\ddot{s}_i = f_i(t, s_1, \ldots, s_4, \dot{s}_1, \ldots, \dot{s}_4)$ — at least not directly. It is possible to reduce such a system to first order, but that basically amounts to extending the system like I showed above by adding four more state variables to it.

This should be enough for you to figure out exactly how to write your MATLAB code, but I'll just point out one thing that may make it easier. You give a dynamical equation for $\ddot{\mathbf{q}}$, but ultimately that comes from a much simpler equation for $\dot{\Omega}_B$, which is the derivative of the angular velocity in the body frame. That equation and the equation for $\dot{\mathbf{q}}$ can be written very easily in terms of just $\Omega_B$ and $\mathbf{q}$, so you could make those your state variables, instead of $\dot{\mathbf{q}}$ and $\mathbf{q}$. That would eliminate a lot of matrix multiplies and inversions, leading to more efficient and possibly more accurate code.

Answer 2

The earlier answer by Mike may well be what you want, and I don't intend to duplicate that, but I'd like to supplement it with some points that I hope will be helpful.

If you are set on the approach of feeding your system of equations $\ddot{\mathbf{q}}=\ldots$ and $\dot{\mathbf{q}}=\ldots$ into an ODE solver in MATLAB, then I would advise taking care and putting in some checks that you've got the correct form of the equations.

Firstly, $\mathbf{q}$ should preserve its normalization, $q_0^2+q_1^2+q_2^2+q_3^2=1$. It's worth checking, as the solver moves forward in time. This also means that $\dot{q}_0 q_0+\dot{q}_1q_1+\dot{q}_2q_2+\dot{q}_3q_3=0$. You need to ensure this is true when you specify your initial conditions. If you calculate $\dot{\mathbf{q}}$ from the kinematics equation, this should be guaranteed automatically, but it does no harm to check.

Secondly, in the absence of an external torque, the rotational kinetic energy and the angular momentum vector $\mathbf{L}=\mathbf{J}\boldsymbol{\omega}$ (calculated in the space-fixed frame) should both be conserved. Again, it's worth checking.

Lastly, I'd like to support the point that Mike made about adopting a simpler approach, based around the body-fixed angular velocity. There are plenty of algorithms for integrating the rotational equations of motion. An approach called "operator splitting" is quite simple to explain, exactly conserves angular momentum, and is symplectic (which is good for energy conservation). Basically, the rotation experienced by the rigid body during a timestep $\delta t$ is split into successive rotations (of duration $\delta t$ or $\frac{1}{2}\delta t$) about each of the body-fixed principal axes (the axis system in which $\mathbf{J}$ is diagonal). A particular pattern of rotations is chosen to ensure time reversibility. The amount of each rotation is determined by the time duration of that step, and the component of angular velocity about the corresponding body-fixed axis. This in turn is determined by the scalar product of $\mathbf{L}$ with that axis vector, divided by the appropriate inertia component. Actually implementing the rotations is fairly simple, using quaternions. Details may be found in Simulating Hamiltonian Mechanics by B Leimkuhler and S Reich, and in papers such as Dullweber, Leimkuhler and McLachlan, J Chem Phys 107, 5840 (1997). The formalism behind these methods may look daunting, but the algorithm is much more straightforward than you might think at first sight.