R(φ, n). Thus four parameters φ, nx , ny , nz describe a rotation, with the

constraint that |n| = 1. Note that R(’φ, ’n) and R(φ, n) represent the

same rotation. The Euler“Rodrigues parameters are expressions of these

four parameters:

[q] = [cos 1 φ, n sin 1 φ]. (15.79)

2 2

They are indeed quaternions that obey the multiplication rule (15.78). They

should be viewed as rotation operators, the product [a][b] meaning the se-

quential operation of ¬rst [b], then [a]. Beware that [q] rotates a vector in a

given coordinate system and not the coordinate system itself, which implies

that [q] is to be identi¬ed with RT and not with R. The unit quaternions

now have the meaning of the identity operator ([1, 0, 0, 0]), and rotations by

π about the x-, y- and z-axes, respectively. Such unit rotations are called

binary rotations. Note that the equivalent rotations (’φ, ’n) and (φ, n) are

given by the same quaternion. Note also that the full range of all normal-

ized quaternions (’1 < qi ¤ +1) includes all possible rotations twice. The

inverse of a normalized quaternion obviously is a rotation about the same

axis but in opposite direction:

[q0 , Q]’1 = [q0 , ’Q]. (15.80)

14 We shall use the notation [q] or [q0 , q1 , q2 , q3 ] or [q,Q] for quaternions.

15 Euler™s theorem: “Two arbitrarily oriented orthonormal bases with common origin P can be

made to coincide with one another by rotating one of them through a certain angle about an

axis which is passing through P and which has the direction of the eigenvector n of the rotation

matrix” (Wittenburg, 1977).

15.7 Rigid body motion 415

r´

r´ φ

φ

φ /2

r

u

nx(nx r)

nxr

n r φ/2

Figure 15.3 Rotation of the vector r to r by a rotation over an angle φ about an

axis n.

We now seek the relation between quaternions and the rotation matrix or

the Euler angles. The latter is rather simple if we realize that the rotation

expressed in Euler angles (•, ‘, ψ) is given by

R = R(ψ, k)R(‘, i)R(•, k) (15.81)

= [cos 1 ψ, 0, 0, sin 1 ψ] [cos 1 ‘, sin 1 ‘, 0, 0] [cos 1 •, 0, 0, sin 1 •]

2 2 2 2 2 2

= [q0 , q1 , q2 , q3 ] (15.82)

with

q0 = cos 1 ‘ cos 1 (• + ψ), (15.83)

2 2

q1 = sin 1 ‘ cos 1 (• ’ ψ), (15.84)

2 2

q2 = sin 1 ‘ cos 1 (• ’ ψ), (15.85)

2 2

q3 = cos 1 ‘ cos 1 (• + ψ). (15.86)

2 2

The relation with the rotation matrix itself can be read from the transfor-

mation of an arbitrary vector r to r by a general rotation over an angle φ

about an axis n, given by the quaternion (15.79). We must be careful with

the de¬nition of rotation: here we consider a vector to be rotated in a ¬xed

coordinate system, while R de¬nes a transformation of the components of

a ¬xed vector to a rotated system of coordinate axes. Each rotation is the

transpose of the other, and the present vector rotation is given by r = RT r.

Refer to Fig. 15.3. De¬ne a vector u, which can be constructed from the

two perpendicular vectors n — (n — r) and n — r:

u = sin 1 φ (n — (n — r)) sin 1 φ + (n — r) cos 1 φ .

2 2 2

416 Lagrangian and Hamiltonian mechanics

With

r = r + 2u = RT r

the rotation matrix can be written out and we obtain

⎛2 ⎞

q0 + q 1 ’ q 2 ’ q 3 2(q1 q3 ’ q0 q2 )

2 2 2 2(q1 q2 + q0 q3 )

R = ⎝ 2(q1 q2 ’ q0 q3 ) 2(q2 q3 + q0 q1 ) ⎠ .

q0 ’ q1 + q2 ’ q3

2 2 2 2

2(q3 q2 ’ q0 q1 ) q0 ’ q1 ’ q2 + q3

2 2 2 2

2(q3 q1 + q0 q2 )

(15.87)

Note Two remarks will be made. The ¬rst concerns the reverse relation, from

R to [q], and the second concerns the symmetry of [q] and its use in generating

random rotations.

The reverse relation is determined by the fact that R has the eigenvalues 1,

exp(iφ) and exp(’iφ), and that the eigenvector belonging to the eigenvalue 1 is n.

So, once these have been determined, [q] follows from (15.79).

As can be seen from (15.87), the four quantities in [q] play an almost, but not

quite symmetric role: q0 di¬ers from q1 , q2 , q3 . This is important if we would fancy

to generate a random rotation. The correct procedure, using quaternions, would

be to generate a vector n randomly distributed over the unit sphere, then generate

a random angle φ in the range (0, 2π), and construct the quaternion from (15.79).

The vector n could be obtained16 by choosing three random numbers x, y, z from

a homogeneous distribution between ’1 and +1, computing r2 = x2 + y 2 + z 2 ,

discarding all triples for which r2 > 1, and scaling all accepted triples to unit

length by division by r. Such a procedure cannot be used for the four dimensions

2

of [q], simply because q0 is di¬erent, e.g., for a random 3D rotation q0 = 1/2

2 2 2

while q1 = q2 = q3 = 1/6. To generate a small random rotation in a given

restricted range, as might be required in a Monte Carlo procedure, ¬rst generate a

random n, then choose a value of φ in the range ”φ as required (both (0, ”φ) and

(’”φ/2, ”φ/2) will do), and construct [q] from (15.79).

What we ¬nally need is the rate of change of [q] due to the angular velocity

ω. Noting that [q] rotates vectors rather than coordinate axes, so that [q] is

™

equal to RT , which is given by (15.60):

™

RT = RT © ,

™ (15.88)

we can cast this into a quaternion multiplication

[q] = [q] [0, 1 ω ].

™ (15.89)

2

Here we have used the fact that © is the time derivative of an angular

rotation around an axis in the direction of ω; the factor 1/2 comes from the

derivative of sin 1 φ. After working this out (reader, please check!) we ¬nd

2

16 See Allen and Tildesley (1987), Appendix G4.

15.8 Holonomic constraints 417

the time derivative of [q]:

⎛ ⎞ ⎛ ⎞

⎛ ⎞

’q1 ’q2 ’q3

q™0

ωx

⎜ q™1 ⎟ 1 ⎜ q0 ’q3 q2 ⎟

⎜ ⎟ ⎜ ⎟ ⎝ ωy ⎠ .

⎝ q™2 ⎠ = 2 ⎝ q3 (15.90)

q0 ’q1 ⎠

ωz

’q2 q1

q™3 q0