Quaternion evolution

Posted on by arkajad

In this note I will repeat some stuff that I was already explaining in the past, but we will also find some important new stuff.
The time evolution differential equation for the attitude matrix Q(t) is

(1)   \begin{equation*}\dot{Q}=QW(\boldsymbol{\Omega}),\end{equation*}

where \boldsymbol{\Omega}=\boldsymbol{\Omega}(t) is the angular velocity vector in the body frame, and W maps any vector \mathbf{v} into skew-symmetric matrix as follows:

(2)   \begin{equation*}W(\mathbf{v})=\begin{bmatrix}0&-v_3&v_2\\v_3&0&-v_1\\-v_2&v_1&0\end{bmatrix}.\end{equation*}

It has the property that for any vector \mathbf{r} we have

(3)   \begin{equation*} W(\mathbf{v})\mathbf{r}=\mathbf{v}\times\mathbf{r}.\end{equation*}

Equation (1) follows immediately from the definition of the angular velocity. The angular velocity \boldsymbol{\omega} in the laboratory frame is defined by

(4)   \begin{equation*}W(\boldsymbol{\omega})=\dot{Q}Q^*.\end{equation*}

But Q maps coordinates of vectors in the body frame into their coordinates in the laboratory frame. Thus \boldsymbol{\omega}=Q\boldsymbol{\Omega}, and therefore

(5)   \begin{equation*}W(Q\boldsymbol{\Omega})=\dot{Q}Q^*.\end{equation*}

But from Eq. (3) we have W(Q\boldsymbol{\Omega})=QW(\boldsymbol{\Omega}) Q^*, therefore

    \[QW(\boldsymbol{\Omega}) Q^*=\dot{Q}Q^*,\]

and so Eq. (1) follows.
While description of rotations in terms of orthogonal 3\times3 matrices in principle suffices in classical mechanics, sometimes it is convenient to use the group of quaternions of unit norm, the group isomorphic to the group \mathrm{SU}(2) used in quantum mechanical description of half-integer spin particles. Quaternions have some advantages in numerical procedures (as for instance in 3D computer games), but they are also convenient for graphical representation of the intrinsic geometry of the rotation group. And this is what interests us, when we plot trajectories representing history of a spinning rigid body.


Quaternions of unit norm, representing rotations, form a 3-dimensional sphere in 4-dimensional Euclidean space, and we are projecting stereographically this sphere onto our familiar three-dimensional space, where we orient ourselves in a usual way known from everyday experience.
The question therefore arises: how the equation describing the time evolution looks like when represented in the quaternion setting?
To derive it we have to return to the fundamental relation between quaternions and rotations of vectors in space.
For every vector \mathbf{v}, with components (v_1,v_2,v_3) denote by \hat{\mathbf{v}} the pure imaginary quaternion defined as:

(6)   \begin{equation*}\hat{\mathbf{v}}=v_1\mathbf{i}+v_2\mathbf{j}+v_3\mathbf{k}.\end{equation*}

Then to unit quaternion q, that is such that qq^*=1, there corresponds rotation matrix R=R(q) such that for all \mathbf{v} the following identity holds

(7)   \begin{equation*}q\hat{\mathbf{v}}q^*=\widehat{R\mathbf{v}}.\end{equation*}

The second interesting property of the hat map is that for all \mathbf{v},\mathbf{w} we have

(8)   \begin{equation*}[\hat{\mathbf{v}},\hat{\mathbf{w}}]=2\widehat{\mathbf{v}\times\mathbf{w}},\end{equation*}

where [\cdot,\cdot] is the commutator. The property follows directly form the definitions and from the quaternion multiplication rules. Every pure imaginary quaternion is of the form \hat{\mathbf{v}} for some \mathbf{v}.

Derivation:
As shown on plaque by the Royal Canal at Broome Bridge in Dublin,

the essence of the quaternions is contained in the four equations

(9)   \begin{equation*}\mathbf{i}^2=\mathbf{j}^2=\mathbf{k}^2=-1,\quad \mathbf{i}\mathbf{j}\mathbf{k}=-1.\end{equation*}

From these, multiplying from the left and/or from the right by \mathbf{i} or \mathbf{j} or \mathbf{k} we derive

(10)   \begin{eqnarray*} \mathbf{i}\mathbf{j}&=&\mathbf{k},\quad \mathbf{j}\mathbf{i}=-\mathbf{k},\notag\\ \mathbf{j}\mathbf{k}&=&\mathbf{i},\quad \mathbf{k}\mathbf{j}=-\mathbf{i},\\ \mathbf{k}\mathbf{i}&=&\mathbf{j},\quad \mathbf{i}\mathbf{k}=-\mathbf{j}\notag. \end{eqnarray*}

From the definition of the hat map we have

(11)   \begin{eqnarray*}\hat{\mathbf{v}}&=&v_1\mathbf{i}+v_2\mathbf{j}+v_3\mathbf{k},\\ \hat{\mathbf{w}}=w_1\mathbf{i}+w_2\mathbf{j}+w_3\mathbf{k}. \end{eqnarray*}

Therefore

    \[ [\hat{\mathbf{v}},\hat{\mathbf{w}}]=\hat{\mathbf{v}}\hat{\mathbf{w}}-\hat{\mathbf{w}},\hat{\mathbf{v}}\]

    \[=(v_1\mathbf{i}+v_2\mathbf{j}+v_3\mathbf{k})(w_1\mathbf{i}+w_2\mathbf{j}+w_3\mathbf{k})-(w_1\mathbf{i}+w_2\mathbf{j}+w_3\mathbf{k})(v_1\mathbf{i}+v_2\mathbf{j}+v_3\mathbf{k}).\]

When we expand the parenthesis and use the multiplication rules in Eqs. (9),(10), we find that the terms with squares cancel out, while the terms with products like \mathbf{i}\mathbf{j} etc. can be organized as follows

    \[ [\hat{\mathbf{v}},\hat{\mathbf{w}}]=2\left[(v_2w_3-w_2v_3)\mathbf{i}+(v_3w_1-w_3v_1)\mathbf{j}+(v_1w_2-w_2v_1)\mathbf{k}\right].\]

The coefficients in the parentheses on the right are now exactly the components of the cross product \mathbf{v}\times\mathbf{w}.
Thus

    \[[\hat{\mathbf{v}},\hat{\mathbf{w}}]=2\widehat{\mathbf{v}\times\mathbf{w}}.\]

Suppose now that q(t) is a trajectory in the space of unit quaternions, while Q(t) is the corresponding trajectory in the rotation group. Thus we have

(12)   \begin{equation*}q(t)\hat{\mathbf{v}}q(t)^*=\widehat{Q(t)\mathbf{v}},\end{equation*}

for all t and all \mathbf{v}. We now differentiate both sides with respect to t. On the left we use the standard product rule for differentiation, but we pay attention so as to preserve the order, because multiplication of quaternions is non commutative. On the right we enter with the differentiation under the “hat”, because it is a linear operation. We obtain:

(13)   \begin{equation*}\dot{q}\hat{\mathbf{v}}q^*+q\hat{\mathbf{v}}\dot{q}^*=\widehat{\dot{Q}\mathbf{v}}.\end{equation*}

Assume now that Q(t) is a solution of Eq. (1), so that \dot{Q}=QW(\boldsymbol{\Omega}). We have

(14)   \begin{equation*}\dot{q}\hat{\mathbf{v}}q^*+q\hat{\mathbf{v}}\dot{q}^*=\widehat{(QW(\boldsymbol{\Omega})\mathbf{v})}.\end{equation*}

We now use Eq. (12) on the right to obtain

(15)   \begin{equation*}\dot{q}\hat{\mathbf{v}}q^*+q\hat{\mathbf{v}}\dot{q}^*=q\widehat{(W(\boldsymbol{\Omega})\mathbf{v})}q^*.\end{equation*}

In order to get rid of \dot{q}^* we differentiate the defining equation of unit quaternions qq^*=1 to obtain

(16)   \begin{equation*}\dot{q}q^*+q\dot{q}^*=0,\end{equation*}

or

(17)   \begin{equation*}\dot{q}^*=-q^*\dot{q}q^*.\end{equation*}

Putting this into Eq. (15) and multiplying both sides on the right by q we obtain

(18)   \begin{equation*}\dot{q}\hat{\mathbf{v}}-q\hat{\mathbf{v}}q^*\dot{q}=q\widehat{(W(\boldsymbol{\Omega})\mathbf{v})}.\end{equation*}

Multiplying from the left by q^*:

(19)   \begin{equation*}q^*\dot{q}\hat{\mathbf{v}}-\hat{\mathbf{v}}q^*\dot{q}=\widehat{(W(\boldsymbol{\Omega})\mathbf{v})}.\end{equation*}

From Eq. (16) it follows that the quaternion q^*\dot{q} is pure imaginary: (q^*\dot{q})^*=-q^*\dot{q}. Therefore there exists vector \mathbf{w} such that

(20)   \begin{equation*}q^*\dot{q}=\hat{\mathbf{w}}.\end{equation*}

Eq. (19) can now be written as

(21)   \begin{equation*}[\hat{\mathbf{w}},\hat{\mathbf{v}}]=\widehat{\boldsymbol{\Omega}\times\mathbf{v}},\end{equation*}

where on the right we have used Eq. (3). Applying Eq. (8) to the left we end with

(22)   \begin{equation*}2\widehat{\mathbf{w}\times\mathbf{v}}=\widehat{\boldsymbol{\Omega}\times\mathbf{v}}.\end{equation*}

Since the above holds for any \mathbf{v}, we deduce that

(23)   \begin{equation*}\hat{\mathbf{w}}=\frac{1}{2}\hat{\boldsymbol{\Omega}}.\end{equation*}

Using (20) we finally obtain:

Evolution equation for quaternions

(24)   \begin{equation*} \dot{q}=\frac{1}{2}q\hat{\boldsymbol{\Omega}}.\end{equation*}

Note: It may be that the final formula can be derived in a shorter way, but I do not know how. It is this last formula that I was using when verifying that the algorithm provided in Meeting with remarkable circles gives indeed a solution of the evolution equation.

13 thoughts on “Quaternion evolution”

  4. So as far as I can understand q(t)=q2(t)q1(t) from earlier post.
    I am intrepretering lower case here to be quaternion form and upper case would be Matrix form i.e Q(t)=Q2(t).Q1(t)
    So the final eqn(21) gives 4 choices for d(psi)/dt to solve.( I can’t get latex working properly).
    I picked one. Gives 1 solution. What initial conditions need to be picked. Solution has one constant to solve for?
    I can try posting the work sheet with a dropbox link you wish.

  5. Hi, Can you show how this arises? I don’t see where it comes from.
    The second interesting property of the hat map is that for allmathbf{v},mathbf{w}we have

    begin{equation*}[hat{mathbf{v}},hat{mathbf{w}}]=2widehat{mathbf{v}timesmathbf{w}},end{equation*}

    where[cdot,cdot] is the commutator.

    2. I am getting \frac{t}{2\sqrt{3}} in the arguments. Which is the same as \frac{\sqrt{3}t}{6}. But I am getting plus signs, not minus signs.
      If you can’t find a mistake in your derivation, let me know, and I will prepare a Maple sheet with my derivation.

      1. Oh! I had used negative values of A1 and A3. Comes out the same as yours now.
        The solution I get for psi(t)=\sqrt {3}\ln \left( \tanh \left( 1/12\,\sqrt {3}t \right) +1 \right) -\sqrt {3}\ln \left( \tanh \left( 1/12\,\sqrt {3}t \right) -1 \right) +2\,\arctan \left( 1/3\,\sqrt {3}\tanh \left( 1/12\,\sqrt {3} t \right) \right) +{\it \_C1}
        setting psi(0)=0 doesnt get me close your expression for psi(t)
        What initial condition should I use?

        1. To fix the constant, so that psi(0)=0, integrate the function of s over s, with s=0..t.

          Then simplify the difference between my expression and your expression. You will get:

          (1/2)*t-I*Pi*sqrt(3)+sqrt(3)*ln(tanh((1/12)*t*sqrt(3))-1)-sqrt(3)*ln(tanh((1/12)*t*sqrt(3))+1)

          Using the definition of tanh in terms of exponential, you can then can show by hand that this is zero.

          Try it. If you can’t let me know, and I will show you how.

          Maple, with its simplification algorithms, is not smart enough to simplify it to zero.

Leave a Reply