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}$.

  • $\begingroup$ Do we normally use \mathbf for the hamiltonian and other variables like the quaternions? I've never seen that before; if I recall correctly, it's usually used as a way to show vectors. You should consider changing that to normal math. People frequently use bold text as a substitute for mathjax (though that isn't encouraged), but you seem to be familiar with mathjax. $\endgroup$ – user191954 Jul 29 '18 at 16:14
  • 1
    $\begingroup$ @Chair it is usual (or at least, not unusual) to use bold for quaternions $\mathbf{q}=(q_0,q_1,q_2,q_3)$, similar to vectors. The OP has not made any reference to a hamiltonian; $\mathbf{H}$ here is a matrix. I think that the use of bold here is correct. $\endgroup$ – user197851 Jul 29 '18 at 16:24
  • $\begingroup$ @LonelyProf Strange... I've been reading a lot about quaternions lately and I've never noticed the bold on the quaternion itself. I just checked and the wikipedia page and the pdf I was using for the stuff (commsp.ee.ic.ac.uk/~mandic/Technical_Report/…) use bold only for the unit vectors... I guess the hamiltonian must've been an accident. I saw the tag hamiltonian formalism and made an assumption. A closer look at the definition of $\mathbf{H}$ shows that I was mistaken. But now the angular velocity is also in bold and I'm severely confused. $\endgroup$ – user191954 Jul 29 '18 at 16:30
  • $\begingroup$ Now that I think of it, it doesn't even make sense to put a quaternion in bold: you'd indicate the individual unit vectors for the last 3 components with bold text, but the number itself is analogous to a complex number. You have any literature about that notation? I'm doing some research on quaternions for a school assignment and that'd be a pretty bad thing to mess up. $\endgroup$ – user191954 Jul 29 '18 at 16:34
  • $\begingroup$ I guess that $\boldsymbol{\omega}$ refers to the angular velocity vector, supplemented by a zero in the zeroth position. But I would encourage the OP to edit the question to define the variables explicitly, since notations do vary. $\endgroup$ – user197851 Jul 29 '18 at 16:35

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.

  • 1
    $\begingroup$ Thanks for the detailed answer! As suggested by you and @LonelyProf , I used the body-fixed angular velocity for the the dynamics equation. Further, $\mathbf{\Omega}_B$ along with the state $\mathbf{q}$ was fed into the kinematics equation set that worked like a charm! The use of $\mathbf{\ddot{q}}$ was motivated by my plan to subsequently incorporate a control strategy that is dependent on the inertial states. $\endgroup$ – Abin Alex Pothen Jul 31 '18 at 14:39

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.

  • $\begingroup$ Regarding the normalization, I actually wrote a paper on integration of angular velocity (the kinematical equation in the question), and found that it's better to relax the constraint. If you rotate vectors as $\mathbf{v}' = \mathbf{q}\, \mathbf{v}\, \mathbf{q}^{-1}$ (using the inverse instead of just the conjugate), then the quaternion just has to have nonzero norm. It turns out that the evolution equation is the same, but inevitable numerical errors in the norm don't impact the results. $\endgroup$ – Mike Jul 31 '18 at 13:19
  • $\begingroup$ @Mike that's interesting; I hadn't thought of it that way. I am still more comfortable, though, keeping $\mathbf{q}$ normalized, and using the conjugate of $\mathbf{q}$ in that formula, or equivalently using the "standard" formula for the rotation matrix (quadratic in the $q_i$). And for the algorithm I mentioned, each step is a simple rotation, so any normalization drift really is at the level of numerical precision rather than algorithm discretization error. For the OP's question, I mainly wanted to suggest checks against mistakes, or unsuitable black-box integration of the ODEs $\endgroup$ – user197851 Jul 31 '18 at 14:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.