Introduction #

... quaternions appear to exude an air of nineteenth century decay, as a rather unsuccessful species in the struggle-for-life of mathematical ideas. Mathematicians, admittedly, still keep a warm place in their hearts for the remarkable algebraic properties of quaternions but, alas, such enthusiasm means little to the harder-headed physical scientist.

— Simon L. Altmann, Rotations, Quaternions, and Double Groups

This article is a work in progress. I'll and refine it a few times over the next week. But until then I don't guarantee correctness.

Early in physics, we learn to model the motion of point particles. We do this by assigning Euclidean coordinates x,y,zx,y,z, and finding differential equations for their evolution. However, in everyday life we rarely encounter point particles. Instead we are surrounded by rigid bodies. While these have a centre of mass that behaves like a point particle, they can also do things like twist and spin. To model these, we need coordinates that can describe all of their degrees of freedom.

In this article we will derive the Euler-Lagrange equations for a rigid body in terms of the quarternions. Euler's equations were re-written in quaternion form in [CR04]. Here we will use [LM18] to generalise this to Lagrangian mechanics on the quaternions.

All images were created by me in Blender. You can use them for your own purposes, but please reference this article. And if you want higher-resolution images feel free to email me.

Euler angles, and their failures #

As before, we will have centre of mass coordinates x,y,zx,y,z. On top of these, we need some more coordinates to describe the orientation of the body. One way of doing this is to use Euler angles. These are three angles representing a series of rotations, that transform the body from its initial orientation to its current one: [B08]

In the image above, θ,ϕ,ψ\theta,\phi,\psi are three angles which transform the coordinate axies x,y,zx,y,z into axies X,Y,ZX,Y,Z representing the current orientation of the body.

I won't go into it here, but we can find differential equations for all of the coordinates x,y,z,θ,ϕ,ψx,y,z,\theta,\phi,\psi , which allow us to solve for how the rigid body evolves in time. These are quite complicated. But even worse, they are singular. If we try and numerically simulate them, there are particular values for the Euler angles where if they are ever reached, the simulation will blow up. This is because, while θ,ϕ,ψ\theta,\phi,\psi should describe a three-dimensional parameter space, there are particular configurations where it collapses to two dimensions. In this situation the coordinates can no longer describe the full range of motion of the object, and the numerics do not survive.

In fact, there is an even more fundamental reason for why they Euler angles, which comes from geometry. The orientation of a rigid body has a naturally spherical geometry. When we write down three angles θ,ϕ,ψ\theta,\phi,\psi, we are using Euclidean flat space to parameterise this curved geometry. Such a thing is an act of violence against the coordinate system, and will always result in a singularity somewhere.

If we want to avoid singularities, we need to find a way to parameterise rigid bodies using a space whose geometry is also curved. We'll see that that the unit quaternions are up to the task.

Crash course in quaternions #

With complex numbers, we add a number ii with the rule i2=−1i^2=-1. With the quaternions, we add three numbers i,j,ki,j,k. These also square to zero:
i2=j2=k2=−1.i^2=j^2=k^2=-1.
Moreover, their products with each other are
ij=k,  jk=i,  ki=j.ij=k,\;jk=i,\;ki=j.
These anti-commute with each other:
ji=−k,  kj=−i,  ik=−j.ji=-k,\;kj=-i,\;ik=-j.
These rules are all you need to know to do computations with the quaternions.

The conjugate of a quaternion is defined similarly to the complex numbers. If
q=q0+iq1+jq2+kq3,q=q_0+iq_1+jq_2+kq_3,
then
q∗=q0−iq1−jq2−kq3.q^*=q_0-iq_1-jq_2-kq_3.
The norm of the quaternion is then
∥q∥2=q∗q=q02+q12+q22+q32.\lVert q\rVert^2=q^*q=q_0^2+q_1^2+q_2^2+q_3^2.
If qq is a unit quaternion, then q−1=q∗.q^{-1}=q^*.

We'll now give a few identities without proof.

Mutliplication rule #

[eqQuaternionMultiplication]: (q1,q⃗1)(q2,q⃗2)=(q1q2−q⃗1⋅q⃗2,  q1q⃗2+q2q⃗1+q⃗1×q⃗2).(q_1,\vec{q}_1)(q_2,\vec{q}_2)=\left(q_1q_2-\vec{q}_1\cdot\vec{q}_2,\;q_1\vec{q}_2+q_2\vec{q}_1+\vec{q}_1\times\vec{q}_2\right).
To see this just expand everything out.

Cross product #

If u⃗,v⃗\vec{u},\vec{v} are pure quaternions, then
[eqCrossProduct]: u⃗v⃗−v⃗u⃗=2u⃗×v⃗.\vec{u}\vec{v}-\vec{v}\vec{u}=2\vec{u}\times\vec{v}.
This follows from [eqQuaternionMultiplication]

Rotations #

Secondly, suppose u^\hat{u} is a unit vector. Then the quaternion
q=(cos⁡(θ/2),sin⁡(θ/2)u^)q=\left(\cos(\theta/2),\sin(\theta/2)\hat{u}\right)
represents a rotation about the axis u^\hat{u} by an angle θ\theta. What this means is that if v⃗\vec{v} is a pure quaternion, and
w⃗=qv⃗q−1,\vec{w}=q\vec{v}q^{-1},
then w⃗\vec{w} is a pure quaternion corresponding to the rotated v⃗\vec{v}. For proof see section 5 of [CR04].

Inner product #

The inner product on the quaternions is defined as
⟨p,q⟩=12(p∗q+q∗p).\langle p,q\rangle=\frac{1}{2}\left(p^*q+q^*p\right).
Suppose p=(p0,p⃗)p=(p_0,\vec{p}), q=(q0,q⃗)q=(q_0,\vec{q}). Then
⟨p,q⟩=Re(pq),\langle p,q\rangle=\mathrm{Re}(pq),
=p0q0+p⃗⋅q⃗.=p_0q_0+\vec{p}\cdot\vec{q}.

And
⟨p,uq⟩=12(p∗uq+q∗u∗p),\langle p,uq\rangle=\frac{1}{2}\left(p^*uq+q^*u^*p\right),
=12((u∗p)∗q+q∗(u∗p)),=\frac{1}{2}\left((u^*p)^*q+q^*(u^*p)\right),
=⟨u∗p,q⟩.=\langle u^*p,q\rangle.

Quaternion description of a rigid body #

We will represent the orientation of the rigid body with a unit quaternion:
q=(cos⁡(θ/2),sin⁡(θ/2)u^).q=\left(\cos(\theta/2),\sin(\theta/2)\hat{u}\right).
Note that we have a one-to-one correspondence between orientations and unit quaternions.

As well as the position, we need the velocity. This is a tangent vector at qq. Let q(t)q(t) be a curve on the sphere, then
q(t)∗q(t)=1.q(t)^*q(t)=1.
Then differentiating both sides, and labelling v=q˙v=\dot{q} as the velocity, gives
[eqTangentCondition]: v∗q+q∗v=0.v^*q+q^*v=0.
Any quaternion vv satisfying this equation is a tangent vector at qq, representing some angular velocity relative to the configuration qq.

Let's try and interpret vv physically, following section 8.2 of [CR04]. Suppose
q=q0+q⃗.q=q_0+\vec{q}.
Then
q∗=q0−q⃗,q^*=q_0-\vec{q},
and
v=v0+v⃗,v=v_0+\vec{v},
where v0,v⃗v_0,\vec{v} are the derivatives of q0,q⃗q_0,\vec{q} with respect to time. The scalar part of q∗vq^*v is
q0v⃗0+q⃗0⋅v⃗0=q0dq0dt+q⃗⋅dq⃗dt,q_0\vec{v}_0+\vec{q}_0\cdot\vec{v}_0=q_0\frac{\mathrm{d}q_0}{\mathrm{d}t}+\vec{q}\cdot\frac{\mathrm{d}\vec{q}}{\mathrm{d}t},
=12ddt∥q∥2,=\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\lVert q\rVert^2,
=0,=0,
where the last line follows since qq is a unit quaternion with constant norm. Thus q∗vq^*v, and its conjugate v∗qv^*q, are pure quaternions. Let's define this pure quaternion as
ω=2q∗v.\omega=2q^*v.
Note that the factor of 1/21/2 is added for convenience based on what will come later. Then
v=12ωq.v=\frac{1}{2}\omega q.

Position xx in the lab frame relates to the body frame XX as
x=qXq∗.x=qXq^*.
Then
x˙=vXq∗+qXv∗,\dot{x}=vXq^*+qXv^*,
=12(ωqXq∗+qXq∗ω∗),=\frac{1}{2}\left(\omega qXq^*+qXq^*\omega^*\right),
=12(ωx−xω),=\frac{1}{2}\left(\omega x-x\omega\right),
where the last line follows since ω\omega is a pure quaternion. From [eqCrossProduct] we find
x˙=ω⃗×x⃗.\dot{x}=\vec{\omega}\times\vec{x}.
Thus ω⃗\vec{\omega} is the angular velocity of the rigid body.

Thus the velocity relates to position as
[eqQuaternionEOM]: q˙=12ωq.\dot{q}=\frac{1}{2}\omega q.

Euler-Lagrange equations #

Suppose we write down a Lagrangian in terms of the quaternion position qq and angular velocity ω\omega. What are the Euler-Lagrange equations, that give us the equations of motion? This will require us to develop variational calculus on the quaternions. We will do this by applying the Lie group formalism in section 8.6 of [LM18].

Infinitesimal variations #

Suppose q(t)q(t) is a path, and let qϵq^{\epsilon} be a variation of this. We need a way of parameterising this in terms of a Euclidean variation. Looking at [eqQuaternionEOM] we can write
qϵ(t)=q(t)eϵη(t),q^{\epsilon}(t)=q(t)e^{\epsilon\eta(t)},
where η(t)\eta(t) is an arbitrary pure quaternion such that η(t0)=η(tf)=0\eta(t_0)=\eta(t_f)=0. The corresponding tangent vector is
q˙(ϵ)(t)=ddtqϵ(t),\dot{q}^{(\epsilon)}(t)=\frac{\mathrm{d}}{\mathrm{d}t}q^{\epsilon}(t),
=q˙(t)eϵη(t)+ϵq(t)η˙(t)eϵη(t).=\dot{q}(t)e^{\epsilon\eta(t)}+\epsilon q(t)\dot{\eta}(t)e^{\epsilon\eta(t)}.
The angular velocity is then
ωϵ(t)=2q∗ϵq˙ϵ,\omega^{\epsilon}(t)=2q^{*\epsilon}\dot{q}^{\epsilon},
=2e−ϵη(t)q∗(t)q˙(t)eϵη(t)+2ϵe−ϵη(t)η˙(t)eϵη(t),=2e^{-\epsilon\eta(t)}q^{*}(t)\dot{q}(t)e^{\epsilon\eta(t)}+2\epsilon e^{-\epsilon\eta(t)}\dot{\eta}(t)e^{\epsilon\eta(t)},
where we have used η(t)∗=−η(t)\eta(t)^*=-\eta(t) since η(t)\eta(t) is a pure quaternion.
This gives us
ωϵ(t)=e−ϵη(t)ω(t)eϵη(t)+2ϵe−ϵη(t)η˙(t)eϵη(t).\omega^{\epsilon}(t)=e^{-\epsilon\eta(t)}\omega(t)e^{\epsilon\eta(t)}+2\epsilon e^{-\epsilon\eta(t)}\dot{\eta}(t)e^{\epsilon\eta(t)}.

Now let us calculate the effect of an infinitesimal variation, when ϵ\epsilon is very small. The variation in orientation is
δq=ddϵqϵ∣ϵ=0,\delta q=\left.\frac{\mathrm{d}}{\mathrm{d}\epsilon}q^{\epsilon}\right\rvert_{\epsilon=0},
which gives
[eqDeltaq]: δq=qη.\delta q=q\eta.
The variation in angular momentum is
δω=ddϵωϵ∣ϵ=0,\delta \omega=\left.\frac{\mathrm{d}}{\mathrm{d}\epsilon}\omega^{\epsilon}\right\rvert_{\epsilon=0},
=−η(t)ω(t)+ω(t)η(t)+2η˙.=-\eta(t)\omega(t)+\omega(t)\eta(t)+2\dot{\eta}.
Thus we find
[eqDeltaomega]: δω=[ω,η]+2η˙.\delta\omega=[\omega,\eta]+2\dot{\eta}.

With these, we can parameterise an arbitrary variation using a pure quaternion η\eta. This lets us turn the problem into one of unconstrained optimisation, since η\eta is an arbitrary three-dimensional vector.

Variation in the action #

Now suppose we have a Lagrangian
L(q,ω)=T(q,ω)−U(q).\mathcal{L}(q,\omega)=T(q,\omega)-U(q).
We want to find the trajectory q(t),ω(t)q(t),\omega(t) which minimises the action
S=∫L(q,ω)dt.\mathcal{S}=\int\mathcal{L}(q,\omega)\mathrm{d}t.

Let us consider an arbitrary variation
Sϵ=∫L(qϵ,ωϵ)dt.\mathcal{S}^{\epsilon}=\int\mathcal{L}(q^{\epsilon},\omega^{\epsilon})\mathrm{d}t.
Then
δSϵ=∫[⟨∂L∂q,δq⟩+⟨∂L∂ω,δω⟩]dt.\delta\mathcal{S}^{\epsilon}=\int\left[\left\langle\frac{\partial\mathcal{L}}{\partial q},\delta q\right\rangle+\left\langle\frac{\partial\mathcal{L}}{\partial\omega},\delta\omega\right\rangle\right]\mathrm{d}t.
Into this we substitute [eqDeltaq] and [eqDeltaomega]. The first term is
⟨∂L∂q,δq⟩=⟨∂L∂q,qη⟩,\left\langle\frac{\partial\mathcal{L}}{\partial q},\delta q\right\rangle=\left\langle\frac{\partial\mathcal{L}}{\partial q},q\eta\right\rangle,
=⟨q∗∂L∂q,η⟩.=\left\langle q^*\frac{\partial\mathcal{L}}{\partial q},\eta\right\rangle.
The second term is
⟨∂L∂ω,δω⟩=⟨∂L∂ω,−η(t)ω(t)+ω(t)η(t)+2η˙⟩,\left\langle\frac{\partial\mathcal{L}}{\partial\omega},\delta\omega\right\rangle=\left\langle\frac{\partial\mathcal{L}}{\partial\omega},-\eta(t)\omega(t)+\omega(t)\eta(t)+2\dot{\eta}\right\rangle,
=⟨−∂L∂ωω(t)∗+ω(t)∗∂L∂ω,η(t)⟩+⟨−2ddt(∂L∂ω),η⟩,=\left\langle-\frac{\partial\mathcal{L}}{\partial\omega}\omega(t)^*+\omega(t)^*\frac{\partial\mathcal{L}}{\partial\omega},\eta(t)\right\rangle+\left\langle-2\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial\mathcal{L}}{\partial\omega}\right),\eta\right\rangle,
=⟨[ω∗,∂L∂ω]−2ddt(∂L∂ω),η⟩.=\left\langle\left[\omega^*,\frac{\partial\mathcal{L}}{\partial\omega}\right]-2\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial\mathcal{L}}{\partial\omega}\right),\eta\right\rangle.

Since ω∗=−ω\omega^*=-\omega this gives us the Euler-Lagrange equation for the quaternions:
ddt(∂L∂ω)=12(q∗∂L∂q+[∂L∂ω,ω]).\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial\mathcal{L}}{\partial\omega}\right)=\frac{1}{2}\left(q^*\frac{\partial\mathcal{L}}{\partial q}+\left[\frac{\partial\mathcal{L}}{\partial\omega},\omega\right]\right).

Examples #

Free body #

Let's now use this to derive Euler's equations for a free rigid body. The Lagrangian in Euclidean space is just the kinetic energy
L=12ω⃗TIω⃗,\mathcal{L}=\frac{1}{2}\vec{\omega}^TI\vec{\omega},
where II is the moment of inertial and ω⃗\vec{\omega} is the angular velocity vector. We need to re-write this in terms of quaternions.

Let ω\omega be the pure quaternion corresponding to ω⃗\vec{\omega}. Define II to be the 4×44\times 4 matrix
(000I),\left(\begin{array}{cc} 0 & 0 \\ 0 & \mathcal{I}\end{array}\right),
where I\mathcal{I} is the 3×33\times 3 inertial tensor. Then we have
L=12⟨ω,Iω⟩.\mathcal{L}=\frac{1}{2}\langle\omega,I\omega\rangle.
Writing it out in components gives
∂L∂ω=Iω,\frac{\partial\mathcal{L}}{\partial\omega}=I\omega,
so the left-hand side of Euler's equations is
ddt(∂L∂ω)=Iω˙.\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial\mathcal{L}}{\partial\omega}\right)=I\dot{\omega}.
The right-hand side is
12[∂L∂ω,ω]=∂L∂ω⃗×ω⃗,\frac{1}{2}\left[\frac{\partial\mathcal{L}}{\partial\omega},\omega\right]=\vec{\frac{\partial\mathcal{L}}{\partial\omega}}\times\vec{\omega},
=(Iω⃗)×ω⃗.=(\mathcal{I}\vec{\omega})\times\vec{\omega}.

Thus we get
Iω⃗′=(Iω⃗)×ω⃗.\mathcal{I}\vec{\omega}'=(\mathcal{I}\vec{\omega})\times\vec{\omega}.

Conclusion #

I'm grateful to Ivan Toftul for introducing me to [CR04] and helping me to understand it.

References #

[CR04] E A Coutsias, & L Romero (2004) The quaternions with an application to rigid body dynamics UNM Digital Repository

[LM18] T Lee, M Leok, & N H McClamroch (2018) Global formulations of Lagrangian and Hamiltonian dynamics on manifolds Springer International Publishing