2Rigid Bodies

The polhode rolls without slipping on the herpolhode lying in the invariable plane.

Herbert Goldstein, Classical Mechanics [20], footnote, p. 207.

The motion of rigid bodies presents many surprising phenomena.

Consider the motion of a top. A top is usually thought of as an axisymmetric body, subject to gravity, with a point on the axis of symmetry that is fixed in space. The top is spun and in general executes some complicated motion. We observe that the top usually settles down into an unusual motion in which the axis of the top slowly precesses about the vertical, apparently moving perpendicular to the direction in which gravity is attempting to accelerate it.

Consider the motion of a book thrown into the air.¹ Books have three main axes. If we idealize a book as a brick with rectangular faces, the three axes are the lines through the centers of opposite faces. Try spinning the book about each axis. The motion of the book spun about the longest and the shortest axis is a simple regular rotation, perhaps with a little wobble depending on how carefully it is thrown. The motion of the book spun about the intermediate axis is qualitatively different: however carefully the book is spun about the intermediate axis, it tumbles.

The rotation of the Moon is peculiar in that the Moon always presents the same face to the Earth, indicating that the rotational period and the orbit period are the same. Considering that the orbit of the Moon is constantly changing because of interactions with the Sun and other planets, and therefore its orbital period is constantly undergoing small variations, we might expect that the face of the Moon that we see would slowly change, but it does not. What is special about the face that is presented to us?

A rigid body may be thought of as a large number of constituent particles with rigid constraints among them. Thus the dynamical principles governing the motion of rigid bodies are the same as those governing the motion of any other system of particles with rigid constraints. What is new here is that the number of constituent particles is very large and we need to develop new tools to handle them effectively.

We have found that a Lagrangian for a system with rigid constraints can be written as the difference of the kinetic and potential energies. The kinetic and potential energies are naturally expressed in terms of the positions and velocities of the constituent particles. To write the Lagrangian in terms of the generalized coordinates and velocities we must specify functions that relate the generalized coordinates to the positions of the constituent particles. In the systems with rigid constraints considered up to now these functions were explicitly given for each of the constituent particles and individually included in the derivation of the Lagrangian. For a rigid body, however, there are too many consituent particles to handle each one of them in this way. We need to find means of expressing the kinetic and potential energies of rigid bodies in terms of the generalized coordinates and velocities, without going through the particle-by-particle details.

The strategy is to first rewrite the kinetic and potential energies in terms of quantities that characterize essential aspects of the distribution of mass in the body and the state of motion of the body. Only later do we introduce generalized coordinates. For the kinetic energy, it turns out a small number of parameters completely specify the state of motion and the relevant aspects of the distribution of mass in the body. For the potential energy, we find that for some specific problems the potential energy can be represented with a small number of parameters, but in general we have to make approximations to obtain a representation with a manageable number of parameters.

2.1   Rotational Kinetic Energy

We consider a rigid body to be made up of a large number of constituent particles with mass m_α, position ${\overrightarrow{x}}_{\alpha }$, and velocities ${\dot{\overrightarrow{x}}}_{\alpha }$, with rigid positional constraints among them. The kinetic energy is

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }{\dot{\overrightarrow{x}}}_{\alpha }·{\dot{\overrightarrow{x}}}_{\alpha }.}& \left(2.1\right)\end{array}$$

It turns out that the kinetic energy of a rigid body can be separated into two pieces: a kinetic energy of translation and a kinetic energy of rotation. Let's see how this comes about.

The configuration of a rigid body is fully specified given the location of any point in the body and the orientation of the body. This suggests that it would be useful to decompose the position vectors for the constituent particles as the sum of the vector $\overrightarrow{X}$ to some reference position in the body and the vector ${\overrightarrow{\xi }}_{\alpha }$ from the reference position to the particular constituent element with index α:

$$\begin{array}{cc}{\overrightarrow{x}}_{\alpha }=\overrightarrow{X}+{\overrightarrow{\xi }}_{\alpha }.& \left(2.2\right)\end{array}$$

Along paths, the velocities are related by

$$\begin{array}{cc}{\dot{\overrightarrow{x}}}_{\alpha }=\dot{\overrightarrow{X}}+{\dot{\overrightarrow{\xi }}}_{\alpha }.& \left(2.3\right)\end{array}$$

So in terms of $\dot{\overrightarrow{X}}$ and ${\dot{\overrightarrow{\xi }}}_{\alpha }$ the kinetic energy is

$$\begin{array}{ll}{\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }\left(\dot{\overrightarrow{X}}+{\dot{\overrightarrow{\xi }}}_{\alpha }\right)·\left(\dot{\overrightarrow{X}}+{\dot{\overrightarrow{\xi }}}_{\alpha }\right)}\hfill & \hfill \\ ={\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }}\left(\dot{\overrightarrow{X}}·\dot{\overrightarrow{X}}+2\dot{\overrightarrow{X}}·{\dot{\overrightarrow{\xi }}}_{\alpha }+{\dot{\overrightarrow{\xi }}}_{\alpha }·{\dot{\overrightarrow{\xi }}}_{\alpha }\right).\hfill & \left(2.4\right)\hfill \end{array}$$

If we select the reference position in the body to be its center of mass,

$$\begin{array}{cc}\overrightarrow{X}=\frac{1}{M}{\displaystyle \sum _{\alpha }{m}_{\alpha }{\overrightarrow{x}}_{\alpha },}& \left(2.5\right)\end{array}$$

where $M={\displaystyle {\sum }_{\alpha }{m}_{\alpha }}$ is the total mass of the body, then

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{m}_{\alpha }{\overrightarrow{\xi }}_{\alpha }={\displaystyle \sum _{\alpha }{m}_{\alpha }\left({\overrightarrow{x}}_{\alpha }-\overrightarrow{X}\right)=0.}}& \left(2.6\right)\end{array}$$

So along paths the relative velocities satisfy

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{m}_{\alpha }{\dot{\overrightarrow{\xi }}}_{\alpha }=0.}& \left(2.7\right)\end{array}$$

The kinetic energy is then

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }\dot{\overrightarrow{X}}·\dot{\overrightarrow{X}}+{\displaystyle {\sum }_{\alpha }\frac{1}{2}{m}_{\alpha }{\dot{\overrightarrow{\xi }}}_{\alpha }·{\dot{\overrightarrow{\xi }}}_{\alpha }.}}& \left(2.8\right)\end{array}$$

The kinetic energy is the sum of the kinetic energy of the motion of the total mass at the center of mass

$$\begin{array}{cc}\frac{1}{2}M\dot{\overrightarrow{X}}·\dot{\overrightarrow{X}},& \left(2.9\right)\end{array}$$

and the kinetic energy of rotation about the center of mass

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }{\dot{\overrightarrow{\xi }}}_{\alpha }·{\dot{\overrightarrow{\xi }}}_{\alpha }.}& \left(2.10\right)\end{array}$$

Written in terms of appropriate generalized coordinates, the kinetic energy is a Lagrangian for a free rigid body. If we choose generalized coordinates so that the center of mass position is entirely specified by some of them and the orientation is entirely specified by others, then the Lagrange equations for a free rigid body will decouple into two groups of equations, one concerned with the motion of the center of mass and one concerned with the orientation.

Such a separation might occur in other problems, such as a rigid body moving in a uniform gravitational field, but in general, potential energies cannot be separated as the kinetic energy separates. So the motion of the center of mass and the rotational motion are usually coupled through the potential. Even in these cases, it is usually an advantage to choose generalized coordinates that separately specify the position of the center of mass and the orientation.

2.2   Kinematics of Rotation

The motion of a rigid body about a center of rotation, a reference position that is fixed with respect to the body, is characterized at each moment by a rotation axis and a rate of rotation. Let's elaborate.

We can get from any orientation of a body to any other orientation of the body by a rotation of the body. That this is true is called Euler's theorem on rotations about a point.² We know that rotations have the property that they do not commute: the composition of successive rotations in general depends on the order of operation. Rotating a book about the $\hat{x}$ axis and then about the ẑ axis puts the book in a different orientation than rotating the book about the ẑ axis and then about the $\hat{x}$ axis. Nevertheless, Euler's theorem states that however many rotations have been composed to reach a given orientation, the orientation could have been reached with a single rotation. Try it! We take a book, rotate it this way, then that, and then some other way—then find the rotation that does the job in one step. So a rotation can be specified by an axis of rotation and the angular amount of the rotation.

We can specify the orientation of a body by specifying the rotation that takes the body to its actual orientation from some reference orientation. As the body moves, the rotation that does this changes.

Let q be the coordinate path that we will use to describe the motion of the body. Let $\mathcal{M}(q(t))$ be the rotation that takes the body from the reference orientation to the orientation specified by q(t) (see figure 2.1). Let ${\overrightarrow{\xi }}_{\alpha }\left(t\right)$ be the vector to some constituent particle with the body in the orientation specified by q(t), and let $\overrightarrow{{\xi }_{\alpha }^{\prime }}$ be the vector to the same constituent with the body in the reference orientation. Then

$$\begin{array}{cc}{\overrightarrow{\xi }}_{\alpha }\left(t\right)=\mathcal{M}\left(q\left(t\right)\right)\overrightarrow{{\xi }_{\alpha }^{\prime }}.& \left(2.11\right)\end{array}$$

The constituent vectors $\overrightarrow{{\xi }_{\alpha }^{\prime }}$ do not depend on the configuration, because they are the vectors to the positions of the constituents with the body in a fixed reference orientation.

To compute the kinetic energy we will accumulate the contributions from all of the constituent mass elements. So we need the velocities of the constituents. The positions of the constituent particles, at a given time t, are

$$\begin{array}{cc}{\overrightarrow{\xi }}_{\alpha }\left(t\right)=\mathcal{M}\left(q\left(t\right)\right)\overrightarrow{{\xi }_{\alpha }^{\prime }}=M\left(t\right)\overrightarrow{{\xi }_{\alpha }^{\prime }},& \left(2.12\right)\end{array}$$

where $M=\mathcal{M}\circ q$. The velocity is the time derivative

$$\begin{array}{cc}D{\overrightarrow{\xi }}_{\alpha }\left(t\right)=DM\left(t\right)\overrightarrow{{\xi }_{\alpha }^{\prime }}.& \left(2.13\right)\end{array}$$

Using equation (2.12), we can write

$$\begin{array}{cc}D{\overrightarrow{\xi }}_{\alpha }\left(t\right)=DM\left(t\right){\left(M\left(t\right)\right)}^{-1}{\overrightarrow{\xi }}_{\alpha }\left(t\right).& \left(2.14\right)\end{array}$$

So we have a time-varying linear differential equation that describes the motion of the constituents. Let's look at the multiplier DM(t)(M(t))⁻¹. Since M(t) is a rotation its matrix representation is an orthogonal matrix M(t), with the property ${(\mathbf{M}(t))}^{-1}={(\mathbf{M}(t))}^{\mathcal{T}}$. Because $M(t){(M(t))}^{\mathcal{T}}=I$, its derivative is:

$$\begin{array}{cc}\mathbf{0}=D\left(M{M}^{\mathcal{T}}\right)=DM{M}^{\mathcal{T}}+MD{M}^{\mathcal{T}}.& \left(2.15\right)\end{array}$$

So

$$\begin{array}{cc}DM{M}^{\mathcal{T}}=-{\left(DM{M}^{\mathcal{T}}\right)}^{\mathcal{T}}.& \left(2.16\right)\end{array}$$

We can conclude that $DM{M}^{\mathcal{T}}$ is antisymmetric.

Let u have components (x, y, z). Every 3 × 3 antisymmetric matrix is of the following form:

$$\begin{array}{cc}\mathcal{A}\left(u\right)=\left(\begin{array}{ccc}0& -z& y\\ z& 0& -x\\ -y& x& 0\end{array}\right).& \left(2.17\right)\end{array}$$

Multiplication by this matrix can be interpreted as the operation of cross product with the vector $\overrightarrow{u}$. The vector $\overrightarrow{u}$ has a matrix representation u.

The inverse of the function $\mathcal{A}$ can be applied to any skew-symmetric matrix: we can use ${\mathcal{A}}^{-1}$ to extract the components of u.

We can interpret multiplication by $DM{M}^{\mathcal{T}}$ as a cross product with a vector that we call $\overrightarrow{\omega }$, the angular velocity vector with components ω. So we can write

$$\begin{array}{cc}\mathit{\omega }={\mathcal{A}}^{-1}\left(DM{M}^{\mathcal{T}}\right).& \left(2.18\right)\end{array}$$

In terms of the angular velocity vector, the differential equations for the motion of the constituents (see equation 2.14) are

$$\begin{array}{cc}D{\overrightarrow{\xi }}_{\alpha }\left(t\right)=\overrightarrow{\omega }\left(t\right)\times {\overrightarrow{\xi }}_{\alpha }\left(t\right).& \left(2.19\right)\end{array}$$

If the angular velocity vector for a body is $\overrightarrow{\omega }$ then the velocities of the constituent particles are perpendicular to the vectors to the constituent particles and proportional to the rate of rotation of the body and the distance of the constituent particle from the instantaneous rotation axis:

$$\begin{array}{cc}{\dot{\overrightarrow{\xi }}}_{\alpha }=\overrightarrow{\omega }\times {\overrightarrow{\xi }}_{\alpha }.& \left(2.20\right)\end{array}$$

The components ω′ of the angular velocity vector on the body axes are $\mathit{\omega }\prime =M^{\mathcal{T}}\mathit{\omega }$, so

$$\begin{array}{cc}\mathit{\omega }\prime =M^{\mathcal{T}}{\mathcal{A}}^{-1}\left(DM{M}^{\mathcal{T}}\right).& \left(2.21\right)\end{array}$$

The relationship of the angular velocity vector to the path is a kinematic relationship; it is valid for any path. Thus we can abstract it to obtain the components of the angular velocity at a moment given the configuration and velocity at that moment.

Implementation of angular velocity functions

The following procedure gives the components of the angular velocity as a function of time along the path:

(define (((M-of-q->omega-of-t M-of-q) q) t)
  (define M-on-path (compose M-of-q q))
  (define (omega-cross t)
    (* ((D M-on-path) t)
       (transpose (M-on-path t))))
  (antisymmetric->column-matrix (omega-cross t)))

The procedure omega-cross produces the matrix representation of $\overrightarrow{\omega }\times$. The procedure antisymmetric->column-matrix, which corresponds to the function ${\mathcal{A}}^{-1}$, is used to extract the components of the angular velocity vector from the skew-symmetric $\overrightarrow{\omega }\times$ matrix.

The components of the angular velocity vector on a basis fixed in the body, as a function of time, along the path are

(define (((M-of-q->omega-body-of-t M-of-q) q) t)
  (* (transpose (M-of-q (q t)))
     (((M-of-q->omega-of-t M-of-q) q) t)))

We can get the procedures of local state that give the angular velocity components by abstracting these procedures along arbitrary paths that have given coordinates and velocities. The abstraction of a procedure of a path to a procedure of state is accomplished by Gamma-bar (see section 1.9):

(define (M->omega M-of-q)
  (Gamma-bar
    (M-of-q->omega-of-t M-of-q)))

(define (M->omega-body M-of-q)
  (Gamma-bar
    (M-of-q->omega-body-of-t M-of-q)))

These procedures give the angular velocities as a function of state. We will see them in action after we get some M-of-q's to work with, starting in section 2.7.

2.3   Moments of Inertia

The rotational kinetic energy is the sum of the kinetic energy of each of the constituents of the rigid body. We can rewrite the rotational kinetic energy in terms of the angular velocity vector and certain aggregate quantities determined by the distribution of mass in the rigid body.

Substituting our representation of the relative velocity vectors into the rotational kinetic energy, we obtain

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }{\dot{\overrightarrow{\xi }}}_{\alpha }·{\dot{\overrightarrow{\xi }}}_{\alpha }={\displaystyle {\sum }_{\alpha }\frac{1}{2}{m}_{\alpha }\left(\overrightarrow{\omega }\times {\overrightarrow{\xi }}_{\alpha }\right)·\left(\overrightarrow{\omega }\times {\overrightarrow{\xi }}_{\alpha }\right).}}& \left(2.22\right)\end{array}$$

We introduce an arbitrary spatially fixed rectangular coordinate system with origin at the center of rotation and with basis vectors ê₀, ê₁, and ê₂, with the property that ê₀ × ê₁ = ê₂. The components of $\overrightarrow{\omega }$ on this coordinate system are ω⁰, ω¹, and ω². Rewriting $\overrightarrow{\omega }$ in terms of its components, the rotational kinetic energy becomes

$$\begin{array}{c}\begin{array}{l}{\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }\left(\left({\displaystyle {\sum }_i{\hat{e}}_i{\omega }^i}\right)\times {\overrightarrow{\xi }}_{\alpha }\right)·\left(\left({\displaystyle {\sum }_j{\hat{e}}_j{\omega }^j}\right)\times {\overrightarrow{\xi }}_{\alpha }\right)}\\ =\frac{1}{2}{\displaystyle {\sum }_{ij}{\omega }^i{\omega }^j}{\displaystyle {\sum }_{\alpha }{m}_{\alpha }}\left({\hat{e}}_i\times {\overrightarrow{\xi }}_{\alpha }\right)·\left({\hat{e}}_j\times {\overrightarrow{\xi }}_{\alpha }\right)\\ =\frac{1}{2}{\displaystyle {\sum }_{ij}{\omega }^i{\omega }^j{I}_{ij},}& \left(2.23\right)\end{array}\end{array}$$

with

$$\begin{array}{cc}{I}_{ij}={\displaystyle \sum _{\alpha }{m}_{\alpha }\left({\hat{e}}_i\times {\overrightarrow{\xi }}_{\alpha }\right)·\left({\hat{e}}_j\times {\overrightarrow{\xi }}_{\alpha }\right).}& \left(2.24\right)\end{array}$$

The nine time-dependent quantities I_ij are the components of the inertia tensor with respect to the chosen coordinate system.

Note what a remarkable form the kinetic energy has taken. All we have done is interchange the order of summations, but now the kinetic energy is written as a sum of products of components of the angular velocity vector, which completely specify how the orientation of the body is changing, and the quantity I_ij, which depends solely on the distribution of mass in the body relative to the chosen coordinate system.

We will deduce a number of properties of the inertia tensor. First, we find a somewhat simpler expression for it. The components of the vector ${\overrightarrow{\xi }}_{\alpha }$ are $\left({\xi }_{\alpha }^0,{\xi }_{\alpha }^1,{\xi }_{\alpha }^2\right)$. If we rewrite ${\overrightarrow{\xi }}_{\alpha }$ as a sum over its components and simplify the elementary vector products of basis vectors, we can obtain the components of the inertia tensor. We can arrange the components of the inertia tensor to form the inertia matrix:

$$\begin{array}{cc}I=\left(\begin{array}{ccc}{I}_{00}& I_{01}& I_{02}\\ I_{10}& I_{11}& I_{12}\\ I_{20}& I_{21}& I_{22}\end{array}\right),& \left(2.25\right)\end{array}$$

where

$$\begin{array}{lll}{I}_{00}\hfill & ={\displaystyle \sum _{\alpha }{m}_{\alpha }({({\xi }_{\alpha }^1)}^2+{({\xi }_{\alpha }^2)}^2)}\hfill & \hfill \\ I_{11}\hfill & ={\displaystyle \sum _{\alpha }{m}_{\alpha }({({\xi }_{\alpha }^0)}^2+{({\xi }_{\alpha }^2)}^2)}\hfill & \hfill \\ I_{22}\hfill & ={\displaystyle \sum _{\alpha }{m}_{\alpha }({({\xi }_{\alpha }^0)}^2+{({\xi }_{\alpha }^1)}^2)}\hfill & \hfill \\ I_{ij}\hfill & =-{\displaystyle \sum _{\alpha }{m}_{\alpha }{\xi }_{\alpha }^i{\xi }_{\alpha }^j}\text{for}i\ne j\hfill & (2.26)\hfill \end{array}$$

Note that the inertia tensor has real components and is symmetric: I_jk = I_kj.

We define the moment of inertia about a line by

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{m}_{\alpha }{\left({\xi }_{\alpha }^{\perp }\right)}^2},& \left(2.27\right)\end{array}$$

where ${\xi }_{\alpha }^{\perp }$ is the perpendicular distance from the line to the constituent with index α. The diagonal components of the inertia tensor I_ii are recognized as the moments of inertia about the lines coinciding with the coordinate axes ê_i. The off-diagonal components of the inertia tensor are called products of inertia.

The rotational kinetic energy of a body depends on the distribution of mass of the body solely through the inertia tensor. Remarkably, the inertia tensor involves only second-order moments of the mass distribution with respect to the center of mass. We might have expected the kinetic energy to depend in a complicated way on all the moments of the mass distribution, interwoven in some complicated way with the components of the angular velocity vector, but this is not the case. This fact has a remarkable consequence: for the motion of a free rigid body the detailed shape of the body does not matter. If a book and a banana have the same inertia tensor, that is, the same second-order mass moments, then if they are thrown in the same way the subsequent motion will be the same, however complicated that motion is. The facts that the book has corners and the banana has a stem do not affect the motion except for their contributions to the inertia tensor. In general, the potential energy of an extended body is not so simple and does indeed depend on all moments of the mass distribution, but for the kinetic energy the second moments are all that matter!

Exercise 2.1: Rotational kinetic energy

Show that the rotational kinetic energy can also be written

$$\begin{array}{cc}{T}_{\mathrm{R}}=\frac{1}{2}I{\omega }^2,& \left(2.28\right)\end{array}$$

where I is the moment of inertia about the line through the center of mass with direction $\hat{\omega }$, and ω is the instantaneous rate of rotation.

Exercise 2.2: Steiner's theorem

Let I be the moment of inertia of a body with respect to some given line through the center of mass. Show that the moment of inertia I′ with respect to a second line parallel to the first is

$$\begin{array}{cc}I\prime =I+M{R}^2& \left(2.29\right)\end{array}$$

where M is the total mass of the body and R is the distance between the lines.

Exercise 2.3: Some useful moments of inertia

Show that the moments of inertia of the following objects are as given:

a. The moment of inertia of a sphere of uniform density with mass M and radius R about any line through the center is $\frac{2}{5}M{R}^2$.

b. The moment of inertia of a spherical shell with mass M and radius R about any line through the center is $\frac{2}{3}M{R}^2$.

c. The moment of inertia of a cylinder of uniform density with mass M and radius R about the axis of the cylinder is $\frac{1}{2}M{R}^2$.

d. The moment of inertia of a thin rod of uniform density per unit length with mass M and length L about an axis perpendicular to the rod through the center of mass is $\frac{1}{12}M{L}^2$.

Exercise 2.4: Jupiter

a. The density of a planet increases toward the center. Provide an argument that the moment of inertia of a planet is less than that of a sphere of uniform density of the same mass and radius.

b. The density as a function of radius inside Jupiter is well approximated by

$$\rho \left(r\right)=\frac{M}{R^3}\frac{\sin \left(\pi r/R\right)}{4r/R},$$

where M is the mass and R is the radius of Jupiter. Find the moment of inertia of Jupiter in terms of M and R.

2.4   Inertia Tensor

The representation of the rotational kinetic energy in terms of the inertia tensor was derived with the help of a rectangular coordinate system with basis vectors ê_i. There was nothing special about this particular rectangular basis. So, the kinetic energy must have the same form in any rectangular coordinate system. We can use this fact to derive how the inertia tensor changes if the body or the coordinate system is rotated.

Let's talk a bit about active and passive rotations. The rotation of the vector $\overrightarrow{x}$ by the rotation R produces a new vector $\overrightarrow{x}\prime =R\overrightarrow{x}$. We may write $\overrightarrow{x}$ in terms of its components with respect to some arbitrary rectangular coordinate system with orthonormal basis vectors ê_i: $\overrightarrow{x}$ = x⁰ê₀ + x¹ê₁ + x²ê₂. Let x indicate the column matrix of components x⁰, x¹, and x² of $\overrightarrow{x}$, and R be the matrix representation of R with respect to the same basis. In these terms the rotation can be written x′ = Rx. The rotation matrix R is a real orthogonal matrix.³ A rotation that carries vectors to new vectors is called an active rotation.

Alternatively, we can rotate the coordinate system by rotating the basis vectors, but leave other vectors that might be represented in terms of them unchanged. If a vector is unchanged but the basis vectors are rotated, then the components of the vector on the rotated basis vectors are not the same as the components on the original basis vectors. Denote the rotated basis vectors by ${\hat{e}}_i^{\prime }=R{\hat{e}}_i$. The component of a vector along a basis vector is the dot product of the vector with the basis vector. So the components of the vector $\overrightarrow{x}$ along the rotated basis ${\hat{e}}_i^{\prime }$ are ${(x\prime )}^i=\overrightarrow{x}·{\hat{e}}_i^{\prime }=\overrightarrow{x}·\left(R{\hat{e}}_i\right)=\left(R^{-1}\overrightarrow{x}\right)·{\hat{e}}_i$.⁴ Thus the components with respect to the rotated basis elements are the same as the components of the rotated vector $R^{-1}\overrightarrow{x}$ with respect to the original basis. In terms of components, if the vector $\overrightarrow{x}$ has components x with respect to the original basis vectors ê_i, then the components x′ of the same vector with respect to the rotated basis vectors ${\hat{e}}_i^{\prime }$ are x′ = R⁻¹x, or equivalently x = Rx′. A rotation that actively rotates the basis vectors, leaving other vectors unchanged, is called a passive rotation. For a passive rotation the components of a fixed vector change as if the vector were actively rotated by the inverse rotation.

With respect to the rectangular basis ê_i the rotational kinetic energy is written

$$\begin{array}{cc}\frac{1}{2}{\displaystyle {\sum }_{ij}{\omega }^i{\omega }^j{I}_{ij}.}& \left(2.30\right)\end{array}$$

In terms of matrix representations, the kinetic energy is

$$\begin{array}{cc}\frac{1}{2}{\mathit{\omega }}^{\mathcal{T}}\mathbf{I}\mathit{\omega },& \left(2.31\right)\end{array}$$

where ω is the column of components representing $\overrightarrow{\omega }$.⁵ If we rotate the coordinate system by the passive rotation R about the center of rotation, the new basis vectors are ${\hat{e}}_i^{\prime }=R{\hat{e}}_i$. The components ω′ of the vector $\overrightarrow{\omega }$ with respect to the rotated coordinate system satisfy

$$\begin{array}{cc}\mathit{\omega }=R\mathit{\omega }\prime ,& \left(2.32\right)\end{array}$$

where R is the matrix representation of R. The kinetic energy is

$$\begin{array}{cc}\frac{1}{2}{(\mathit{\omega }\prime )}^{\mathcal{T}}{R}^{\mathcal{T}}IR\mathit{\omega }\prime .& \left(2.33\right)\end{array}$$

However, if we had started with the basis ${\hat{e}}_i^{\prime }$, we would have written the kinetic energy directly as

$$\begin{array}{cc}\frac{1}{2}{(\mathit{\omega }\prime )}^{\mathcal{T}}I\prime \mathit{\omega }\prime ,& \left(2.34\right)\end{array}$$

where the components are taken with respect to the ${\hat{e}}_i^{\prime }$ basis. Comparing the two expressions, we see that

$$\begin{array}{cc}I\prime =R^{\mathcal{T}}IR.& \left(2.35\right)\end{array}$$

Thus the inertia matrix transforms by a similarity transformation.⁶

2.5   Principal Moments of Inertia

We can use the transformation properties of the inertia tensor (2.35) to show that there are special rectangular coordinate systems for which the inertia tensor I′ is diagonal, that is, $I_{ij}^{\prime }=0$ for i ≠ j. Let's assume that I′ is diagonal and solve for the rotation matrix R that does the job. Multiplying both sides of (2.35) on the left by R, we have

$$\begin{array}{cc}RI\prime =IR.& \left(2.36\right)\end{array}$$

We can examine pieces of this matrix equation by multiplying on the right by a trivial column vector that picks out a particular column. So we multiply on the right by the column matrix representation e_i of each of the coordinate unit vectors ê_i. These column matrices have a one in the ith row and zeros otherwise. Using $e_i^{\prime }=R{e}_i$, we find

$$\begin{array}{cc}RI\prime e_i=IR{e}_i=I{e}_i^{\prime }.& \left(2.37\right)\end{array}$$

On the other hand, the matrix I′ is diagonal, so

$$\begin{array}{cc}RI\prime e_i=R{e}_i{I}_{ii}^{\prime }=I_{ii}^{\prime }{e}_i^{\prime }.& \left(2.38\right)\end{array}$$

So, from equations (2.37) and (2.38), we have

$$\begin{array}{cc}I{e}_i^{\prime }=I_{ii}^{\prime }{e}_i^{\prime },& \left(2.39\right)\end{array}$$

which we recognize as an equation for the eigenvalue $I_{ii}^{\prime }$ and $e_i^{\prime }$, the column matrix of components of the associated eigenvector.

From $e_i^{\prime }$ = Re_i, we see that the $e_i^{\prime }$ are the columns of the rotation matrix R. Now, rotation matrices are orthogonal, so $R^{\mathcal{T}}R=1$; thus the columns of the rotation matrix must be orthonormal—that is, ${\left(e_i^{\prime }\right)}^{\mathcal{T}}{e}_j^{\prime }={\delta }_{ij}$, where δ_ij is one if i = j and zero otherwise. But the eigenvectors that are solutions of equation (2.39) are not necessarily even orthogonal. So we are not done yet.

If a matrix is real and symmetric then the eigenvalues are real. Furthermore, if the eigenvalues are distinct then the eigenvectors are orthogonal. However, if the eigenvalues are not distinct then the directions of the eigenvectors for the degenerate eigenvalues are not uniquely determined—we have the freedom to choose particular $e_i^{\prime }$ that are orthogonal.⁷ The linearity of equation (2.39) implies that the $e_i^{\prime }$ can be normalized. Thus whether or not the eigenvalues are distinct we can obtain an orthonormal set of $e_i^{\prime }$. This is enough to reconstruct a rotation matrix R that does the job we asked of it: to rotate the coordinate system to a configuration such that the inertia tensor is diagonal. If the eigenvalues are not distinct, the rotation matrix R is not uniquely defined—there is more than one rotation matrix R that does the job.

The eigenvectors and eigenvalues are determined by the requirement that the inertia tensor be diagonal with respect to the rotated coordinate system. Thus the rotated coordinate system has a special orientation with respect to the body. The basis vectors ${\hat{e}}_i^{\prime }$ therefore actually point along particular directions in the body. We define the axes in the body through the center of mass with these directions to be the principal axes. With respect to the coordinate system defined by ${\hat{e}}_i^{\prime }$, the inertia tensor is diagonal, by construction, with the eigenvalues $I_{ii}^{\prime }$ on the diagonal. Thus the moments of inertia about the principal axes are the eigenvalues $I_{ii}^{\prime }$. We call the moments of inertia about the principal axes the principal moments of inertia.

For convenience, we often label the principal moments of inertia according to their size: A ≤ B ≤ C, with principal axis unit vectors â, $\hat{b}$, ĉ, respectively. The positive direction along the principal axes can be chosen so that â, $\hat{b}$, ĉ form a right-handed rectangular coordinate basis.

Let x represent the matrix of components of a vector $\overrightarrow{x}$ with respect to the basis vectors ê_i. Recall that the components x′ of a vector $\overrightarrow{x}$ with respect to the principal axis unit vectors ${\hat{e}}_i^{\prime }$ satisfy

$$\begin{array}{cc}x\prime =R^{\mathcal{T}}x.& \left(2.40\right)\end{array}$$

The components of a vector on the principal axis basis are sometimes called the body components of the vector.

If we choose the reference orientation of the body so that the principal axes are aligned with the spatial axes $\hat{x}$, ŷ, ẑ, then the rotation R that diagonalizes the inertia matrix becomes the rotation M shown in figure 2.1. The axes â′, $\hat{b}\prime$, ĉ′ then become the principal axes. The rotation matrix M multiplies the column of components of a vector on the principal axes to make a column of components of the vector in space.

Now let's rewrite the kinetic energy in terms of the principal moments of inertia. If we choose our rectangular coordinate system so that it coincides with the principal axes then the calculation is simple. Let the components of the angular velocity vector on the principal axes be (ω^a, ω^b, ω^c). Then, keeping in mind that the inertia tensor is diagonal with respect to the principal axis basis, the kinetic energy is just

$$\begin{array}{cc}{T}_{\mathrm{R}}=\frac{1}{2}\left[A{({\omega }^a)}^2+B{({\omega }^b)}^2+C{({\omega }^c)}^2\right].& (2.41)\end{array}$$

Or as a program:

(define ((T-body A B C) omega-body)
  (* 1/2
     (+ (* A (square (ref omega-body 0)))
        (* B (square (ref omega-body 1)))
        (* C (square (ref omega-body 2))))))

Exercise 2.5: A constraint on the moments of inertia

Show that the sum of any two of the moments of inertia is greater than or equal to the third moment of inertia. You may assume the moments of inertia are with respect to orthogonal axes.

Exercise 2.6: Principal moments of inertia

For each of the configurations described below find the principal moments of inertia with respect to the center of mass, and find the corresponding principal axes.

a. A regular tetrahedron consisting of four equal point masses tied together with rigid massless wire.

b. A cube of uniform density.

c. Five equal point masses rigidly connected by massless stuff. The point masses are at the rectangular coordinates

(−1, 0, 0), (1, 0, 0), (1, 1, 0), (0, 0, 0), (0, 0, 1).

Exercise 2.7: This book

Measure this book. You will admit that it is pretty dense. Don't worry, you will get to throw it later. Show that the principal axes are the lines connecting the centers of opposite faces of the idealized brick approximating the book. Compute the corresponding principal moments of inertia.

2.6   Vector Angular Momentum

The vector angular momentum of a particle is the cross product of its position vector and its linear momentum vector. For a rigid body the vector angular momentum is the sum of the vector angular momentum of each of the constituents. Here we find an expression for the vector angular momentum of a rigid body in terms of the inertia tensor and the angular velocity vector.

The vector angular momentum of a rigid body is

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{\overrightarrow{x}}_{\alpha }\times \left(m_{\alpha }{\dot{\overrightarrow{x}}}_{\alpha }\right),}& \left(2.42\right)\end{array}$$

where ${\overrightarrow{x}}_{\alpha }$, ${\dot{\overrightarrow{x}}}_{\alpha }$, and m_α are the positions, velocities, and masses of the constituent particles. It turns out that the vector angular momentum decomposes into the sum of the angular momentum of the center of mass and the rotational angular momentum about the center of mass, just as the kinetic energy separates into the kinetic energy of the center of mass and the kinetic energy of rotation. As in the kinetic energy demonstration (section 2.1), decompose the position into the vector to the center of mass $\overrightarrow{X}$ and the vectors from the center of mass to the constituent mass elements ${\overrightarrow{\xi }}_{\alpha }$:

$$\begin{array}{cc}{\overrightarrow{x}}_{\alpha }=\overrightarrow{X}+{\overrightarrow{\xi }}_{\alpha },& \left(2.43\right)\end{array}$$

with velocities

$$\begin{array}{cc}{\dot{\overrightarrow{x}}}_{\alpha }=\dot{\overrightarrow{X}}+{\dot{\overrightarrow{\xi }}}_{\alpha }.& \left(2.44\right)\end{array}$$

Substituting, the angular momentum is

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{m}_{\alpha }\left(\overrightarrow{X}+{\overrightarrow{\xi }}_{\alpha }\right)\times \left(\dot{\overrightarrow{X}}+{\dot{\overrightarrow{\xi }}}_{\alpha }\right).}& \left(2.45\right)\end{array}$$

Multiplying out the product, and using the fact that $\overrightarrow{X}$ is the center of mass and $M={\displaystyle {\sum }_{\alpha }{m}_{\alpha }}$ is the total mass of the body, the angular momentum is

$$\begin{array}{cc}\overrightarrow{X}\times \left(M\dot{\overrightarrow{X}}\right)+{\displaystyle \sum _{\alpha }{\overrightarrow{\xi }}_{\alpha }\times \left(m_{\alpha }{\dot{\overrightarrow{\xi }}}_{\alpha }\right).}& \left(2.46\right)\end{array}$$

The angular momentum of the center of mass is

$$\begin{array}{cc}\overrightarrow{X}\times \left(M\dot{\overrightarrow{X}}\right),& \left(2.47\right)\end{array}$$

and the rotational angular momentum is

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{\overrightarrow{\xi }}_{\alpha }\times \left(m_{\alpha }{\dot{\overrightarrow{\xi }}}_{\alpha }\right).}& \left(2.48\right)\end{array}$$

Using ${\dot{\overrightarrow{\xi }}}_{\alpha }=\overrightarrow{\omega }\times {\overrightarrow{\xi }}_{\alpha }$, we get the rotational angular momentum vector

$$\begin{array}{cc}\overrightarrow{L}={\displaystyle \sum _{\alpha }{m}_{\alpha }{\overrightarrow{\xi }}_{\alpha }\times \left(\overrightarrow{\omega }\times {\overrightarrow{\xi }}_{\alpha }\right).}& \left(2.49\right)\end{array}$$

We can also reexpress the rotational angular momentum in terms of the angular velocity vector and the inertia tensor, as we did for the kinetic energy. In terms of components with respect to the basis {ê₀, ê₁, ê₂}, this is

$$\begin{array}{cc}{L}_j={\displaystyle \sum _k{I}_{jk}{\omega }^k},& \left(2.50\right)\end{array}$$

where I_jk are the components of the inertia tensor (2.24). The angular momentum and the kinetic energy are expressed in terms of the same inertia tensor.

With respect to the principal-axis basis, the components of the angular momentum have a particularly simple form:

$$\begin{array}{lll}{L}_a\hfill & =A{\omega }^a\hfill & \hfill \\ L_b\hfill & =B{\omega }^b\hfill & \hfill \\ L_c\hfill & =C{\omega }^c\hfill & \left(2.51\right)\hfill \end{array}$$

Since the angular momenta are the partial derivatives of T_R (see equation 2.41) with respect to the angular velocities, they must be grouped as a down tuple (in matrix language, a row matrix): L′ = [L_a, L_b, L_c]. As a program:

(define ((L-body A B C) omega-body)
  (down (* A (ref omega-body 0))
        (* B (ref omega-body 1))
        (* C (ref omega-body 2))))

If M is the matrix representation of the rotation that takes an angular-velocity vector $\overrightarrow{\omega }\prime$ to a rotated vector $\overrightarrow{\omega }$, the components transform as ω = Mω′.

When working with matrices it is more convenient to work with a column matrix of the angular momentum components, so we introduce $\overline{L}=L^{\mathcal{T}}$. Using ω = Mω′ and equation (2.35) with R replaced by M we derive an expression for the angular momentum

$$\begin{array}{cc}\overline{L}=I\mathit{\omega }=MI\prime \mathit{\omega }\prime =M\overline{L}\prime .& \left(2.52\right)\end{array}$$

Transposing this result, we see that the angular momentum components must transform as $L=L\prime M^{\mathcal{T}}$:

(define (((L-space M) A B C) omega-body)
  (* ((L-body A B C) omega-body)
     (transpose M)))

Exercise 2.8: Rotational angular momentum

Verify that expression (2.50) for the components of the rotational angular momentum (2.49) in terms of the inertia tensor is correct.

2.7   Euler Angles

To go further we must finally specify a set of generalized coordinates. We first do this using the traditional Euler angles. Later, we find other ways of describing the orientation of a rigid body.

We are using an intermediate representation of the orientation in terms of the function $\mathcal{M}$ of the generalized coordinates that gives the rotation that takes the body from some reference orientation and rotates it to the orientation specified by the generalized coordinates. Here we take the reference orientation so that principal-axis unit vectors â, $\hat{b}$, ĉ are coincident with the basis vectors ê_i, labeled here by $\hat{x}$, ŷ, ẑ.

We define the Euler angles in terms of simple rotations about the coordinate axes. Let R_x(ψ) be a right-handed rotation about the $\hat{x}$ axis by the angle ψ, and let R_z(ψ) be a right-handed rotation about the ẑ axis by the angle ψ. The function $\mathcal{M}$ for Euler angles is written as a composition of three of these simple coordinate axis rotations:

$$\begin{array}{cc}\mathcal{M}\left(\theta ,\phi ,\psi \right)=R_z\left(\phi \right)\circ R_x\left(\theta \right)\circ R_z\left(\psi \right),& \left(2.53\right)\end{array}$$

for the Euler angles θ, φ, ψ.

The Euler angles can specify any orientation of the body, but the orientation does not always correspond to a unique set of Euler angles. In particular, if θ = 0 then the orientation is dependent only on the sum φ + ψ, so the orientation does not uniquely determine either φ or ψ.

Exercise 2.9: Euler angles

It is not immediately obvious that all orientations can be represented in terms of the Euler angles. To show that the Euler angles are adequate to represent all orientations, solve for the Euler angles that give an arbitrary rotation R. Keep in mind that some orientations do not correspond to a unique representation in terms of Euler angles.

Though the Euler angles allow us to specify all orientations and thus can be used as generalized coordinates, the definition of Euler angles is pretty arbitrary. In fact no reasoning has led us to them. This is reflected in our presentation of them by just saying “here they are.” Euler angles are well suited for some problems, but cumbersome for others.

There are other ways of defining similar sets of angles. For instance, we could also take our generalized coordinates to satisfy

$$\begin{array}{cc}\mathcal{M}\prime \left(\theta ,\phi ,\psi \right)=R_x\left(\phi \right)\circ R_y\left(\theta \right)\circ R_z\left(\psi \right).& \left(2.54\right)\end{array}$$

Such alternatives to the Euler angles sometimes come in handy.

Each of the fundamental rotations can be represented as a matrix. The rotation matrix representing a right-handed rotation about the ẑ axis by the angle ψ is

$$\begin{array}{cc}{R}_z\left(\psi \right)=\left(\begin{array}{ccc}\cos \psi & -\sin \psi & 0\\ \sin \psi & \cos \psi & 0\\ 0& 0& 1\end{array}\right)& \left(2.55\right)\end{array}$$

and a right-handed rotation about the x axis by the angle ψ is represented by the matrix

$$\begin{array}{cc}{R}_x\left(\psi \right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \psi & -\sin \psi \\ 0& \sin \psi & \cos \psi \end{array}\right).& \left(2.56\right)\end{array}$$

The matrix that represents the rotation that carries the body from its reference orientation to the actual orientation is

$$\begin{array}{cc}M\left(\theta ,\phi ,\psi \right)=R_z\left(\phi \right)R_x\left(\theta \right)R_z\left(\psi \right).& \left(2.57\right)\end{array}$$

The rotation matrices and their product can be constructed by simple programs:

(define (Rz-matrix angle)
  (matrix-by-rows
    (list (cos angle) (- (sin angle)) 0)
    (list (sin angle) (cos angle) 0)
    (list 0 0 1)))

(define (Rx-matrix angle)
  (matrix-by-rows
    (list 1 0 0)
    (list 0 (cos angle) (- (sin angle)))
    (list 0 (sin angle) (cos angle))))

(define (Euler->M angles)
  (let ((theta (ref angles 0))
        (phi (ref angles 1))
        (psi (ref angles 2)))

    (* (Rz-matrix phi)
       (Rx-matrix theta)
       (Rz-matrix psi))))

Now that we have a procedure that implements a sample $\mathcal{M}$, we can find the components of the angular velocity vector and the body components of the angular velocity vector using the procedures M-of-q->omega-of-t and M-of-q->omega-body-of-t from section 2.2. For example,

(show-expression
  (((M-of-q->omega-body-of-t Euler->M)
    (up (literal-function 'theta)
        (literal-function 'phi)
        (literal-function 'psi)))
   't))

To construct the kinetic energy we need the procedure of state that gives the body components of the angular velocity vector:

(show-expression
  ((M->omega-body Euler->M)
   (up 't
       (up 'theta 'phi 'psi)
       (up 'thetadot 'phidot 'psidot))))

We capture this result as a procedure:

(define (Euler-state->omega-body local)
  (let ((q (coordinate local)) (qdot (velocity local)))
    (let ((theta (ref q 0))
          (psi (ref q 2))
          (thetadot (ref qdot 0))
          (phidot (ref qdot 1))
          (psidot (ref qdot 2)))
      (let ((omega-a (+
                       (* thetadot (cos psi))
                       (* phidot (sin theta) (sin psi))))
            (omega-b (+
                       (* -1 thetadot (sin psi))
                       (* phidot (sin theta) (cos psi))))
            (omega-c (+
                       (* phidot (cos theta)) psidot)))
        (up omega-a omega-b omega-c)))))

The kinetic energy can now be written:

(define ((T-body-Euler A B C) local)
  ((T-body A B C)
   (Euler-state->omega-body local)))

We can define procedures to calculate the components of the angular momentum on the principal axes:

(define ((L-body-Euler A B C) local)
  ((L-body A B C)
   (Euler-state->omega-body local)))

We then transform the components of the angular momentum on the principal axes to the components on the fixed basis ê_i:

(define ((L-space-Euler A B C) local)
  (let ((angles (coordinate local)))
    (* ((L-body-Euler A B C) local)
       (transpose (Euler->M angles)))))

These procedures are local state functions, like Lagrangians.

2.8   Motion of a Free Rigid Body

The kinetic energy, expressed in terms of a suitable set of generalized coordinates, is a Lagrangian for a free rigid body. In section 2.1 we found that the kinetic energy of a rigid body can be written as the sum of the rotational kinetic energy and the translational kinetic energy. If we choose one set of coordinates to specify the position and another set to specify the orientation, the Lagrangian becomes a sum of a translational Lagrangian and a rotational Lagrangian. The Lagrange equations for translational motion are not coupled to the Lagrange equations for the rotational motion. For a free rigid body the translational motion is just that of a free particle: uniform motion. Here we concentrate on the rotational motion of the free rigid body. We can adopt the Euler angles as the coordinates that specify the orientation; the rotational kinetic energy was expressed in terms of Euler angles in the previous section.

Conserved quantities

The Lagrangian for a free rigid body has no explicit time dependence, so we can deduce that the energy, which is just the kinetic energy, is conserved by the motion.

The Lagrangian does not depend on the Euler angle φ, so we can deduce that the momentum conjugate to this coordinate is conserved. An explicit expression for the momentum conjugate to φ is

(define Euler-state
  (up 't
      (up 'theta 'phi 'psi)
      (up 'thetadot 'phidot 'psidot)))

(show-expression
  (ref (((partial 2) (T-body-Euler 'A 'B 'C)) Euler-state)
       1))

We know that this complicated quantity is conserved by the motion of the rigid body because of the symmetries of the Lagrangian.

If there are no external torques, then we expect that the vector angular momentum will be conserved. We can verify this using the Lagrangian formulation of the problem. First, we note that L_z is the same as p_φ. We can check this by direct calculation:

(- (ref ((L-space-Euler 'A 'B 'C) Euler-state)
        2)
   (ref (((partial 2) (T-body-Euler 'A 'B 'C)) Euler-state)
        1))

0

We know that p_φ is conserved because the Lagrangian for the free rigid body did not mention φ, so now we know that L_z is conserved. Since the orientation of the coordinate axes is arbitrary, we know that if any rectangular component is conserved then all of them are. So the vector angular momentum is conserved for the free rigid body.

We could have seen this with the help of Noether's theorem (see section 1.8.5). There is a continuous family of rotations that can transform any orientation into any other orientation. The orientation of the coordinate axes we used to define the Euler angles is arbitrary, and the kinetic energy (the Lagrangian) is the same for any choice of coordinate system. Thus the situation meets the requirements of Noether's theorem, which tells us that there is a conserved quantity. In particular, the family of rotations around each coordinate axis gives us conservation of the angular momentum component on that axis. We construct the vector angular momentum by combining these contributions.

Exercise 2.10: Uniformly accelerated rigid body

Show that a rigid body subject to a uniform acceleration rotates as a free rigid body, while the center of mass has a parabolic trajectory.

Exercise 2.11: Conservation of angular momentum

Fill in the details of the argument that Noether's theorem implies that vector angular momentum is conserved by the motion of the free rigid body.

2.8.1 Computing the Motion of Free Rigid Bodies

Lagrange's equations for the motion of a free rigid body in terms of Euler angles are quite disgusting, so we will not show them here. However, we will use the Lagrange equations to explore the motion of the free rigid body.

Before doing this it is worth noting that the equations of motion in Euler angles are singular for some configurations, because for these configurations the Euler angles are not uniquely defined. If we set θ = 0 then an orientation does not correspond to a unique value of φ and ψ; only their sum determines the orientation.

The singularity arises in the explicit Lagrange equations when we attempt to solve for the second derivative of the generalized coordinates in terms of the generalized coordinates and the generalized velocities (see section 1.7). The isolation of the second derivative requires multiplying by the inverse of ∂₂∂₂L. The determinant of this quantity becomes zero when the Euler angle θ is zero:

(show-expression
  (determinant
    (((square (partial 2)) (T-body-Euler 'A 'B 'C))
     Euler-state)))

So when θ is zero, we cannot solve for the second derivatives. When θ is small, the Euler angles can move very rapidly, and thus may be difficult to compute reliably. Of course, the motion of the rigid body is perfectly well behaved for any orientation. This is a problem of the representation of that motion in Euler angles; it is a “coordinate singularity.”

One solution to this problem is to use another set of Euler-like coordinates for which Lagrange's equations have singularities for different orientations, such as those defined in equation (2.54). So if as the calculation proceeds the trajectory comes close to a singularity in one set of coordinates, we can switch coordinate systems and use another set for a while until the trajectory encounters another singularity. This solves the problem, but it is cumbersome. For the moment we will ignore this problem and compute some trajectories, being careful to limit our attention to trajectories that avoid the singularities.

We will compute some trajectories by numerical integration and check our integration process by seeing how well energy and angular momentum are conserved. Then, we will investigate the evolution of the components of angular momentum on the principal axis basis. We will discover that we can learn quite a bit about the qualitative behavior of rigid bodies by combining the information we get from the energy and angular momentum.

To develop a trajectory from initial conditions we integrate the Lagrange equations, as we did in chapter 1. The system derivative is obtained from the Lagrangian:

(define (rigid-sysder A B C)
  (Lagrangian->state-derivative (T-body-Euler A B C)))

The following program monitors the errors in the energy and in the components of the angular momentum:

(define ((monitor-errors win A B C L0 E0) state)
  (let ((t (time state))
        (L ((L-space-Euler A B C) state))
        (E ((T-body-Euler A B C) state)))
    (plot-point win t (relative-error (ref L 0) (ref L0 0)))
    (plot-point win t (relative-error (ref L 1) (ref L0 1)))
    (plot-point win t (relative-error (ref L 2) (ref L0 2)))
    (plot-point win t (relative-error E E0))))

(define (relative-error value reference-value)
  (if (zero? reference-value)
    (error "Zero reference value -- RELATIVE-ERROR")
    (/ (- value reference-value) reference-value)))

We make a plot window to display the errors:

(define win (frame 0.0 100.0 -1.0e-12 1.0e-12))

The default integration method used by the system is Bulirsch–Stoer (bulirsch-stoer), but here we set the integration method to be quality-controlled Runge–Kutta (qcrk4), because the error plot is more interesting:

(set-ode-integration-method! 'qcrk4)

We use evolve to investigate the evolution:

(up 1.0 0.0 0.0)

(up 0.1 0.1 0.1))))

(let ((L0 ((L-space-Euler A B C) state0))

(E0 ((T-body-Euler A B C) state0)))

((evolve rigid-sysder A B C)

state0

(monitor-errors win A B C L0 E0)

The plot that is developed of the relative errors in the components of the angular momenta and the energy (see figure 2.2) shows that we have been successful in controlling the error in the conserved quantities. This should give us some confidence in the trajectory that is evolved.

2.8.2 Qualitative Features of Free Rigid Body Motion

The evolution of the components of the angular momentum on the principal axes has a remarkable property. For almost every initial condition the body components of the angular momentum periodically trace a simple closed curve.

We can see this by investigating a number of trajectories and plotting the components of angular momentum of the body on the principal axes (see figure 2.3). To make this figure a number of trajectories of equal energy were computed. The three-dimensional space of body components is projected onto a two-dimensional plane for display. Points on the back of this projection of the ellipsoid of constant energy are plotted with lower density than points on the front of the ellipsoid. For most initial conditions we find a one-dimensional simple closed curve. Some trajectories on the front side appear to cross trajectories on the back side, but this is an artifact of projection. There is also a family of trajectories that appear to intersect in two points, one on the front side and one on the back side. The curve that is the union of these trajectories is called a separatrix; it separates different types of motion.

What is going on? The state space for a free rigid body is six-dimensional: the three Euler angles and their time derivatives. We know four constants of the motion—the three spatial components of the angular momentum, L_x, L_y, and L_z, and the energy, E. Thus, the motion is restricted to a two-dimensional region of the state space.⁸ Our experiment shows that the components of the angular momentum trace one-dimensional closed curves in the angular-momentum subspace, so there is something more going on here.

The total angular momentum is conserved if all of the components are, so we also have the constant

$$\begin{array}{cc}{L}^2=L_x^2+L_y^2+L_z^2.& \left(2.58\right)\end{array}$$

The spatial components of the angular momentum do not change, but of course the projections of the angular momentum onto the principal axes do change because the axes move as the body moves. However, the magnitude of the angular momentum vector is the same whether it is computed from components on the fixed basis or components on the principal axis basis. So, the combination

$$\begin{array}{cc}{L}^2=L_a^2+L_b^2+L_c^2,& \left(2.59\right)\end{array}$$

is conserved.

Using the expressions (2.51) for the components of the angular momentum in terms of the components of the angular velocity vector on the principal axes, the kinetic energy (2.41) can be rewritten in terms of the angular momentum components on the principal axes:

$$\begin{array}{cc}E=\frac{1}{2}\left(\frac{L_a^2}{A}+\frac{L_b^2}{B}+\frac{L_c^2}{C}\right).& \left(2.60\right)\end{array}$$

The two conserved quantities (2.59 and 2.60) provide constraints on how the components of the angular momentum vector on the principal axes can change. We recognize the conservation of angular momentum constraint (2.59) as the equation of a sphere, and the conservation of kinetic energy constraint (2.60) as the equation for a triaxial ellipsoid. For every trajectory both constraints are satisfied, so the components of the angular momentum move on the intersection of these two surfaces, the energy ellipsoid and the angular momentum sphere. The intersection of an ellipsoid and a sphere with the same center is typically two closed curves, so an orbit is confined to one of these curves. This sheds light on the puzzle posed at the beginning of this section.

Because of our ordering A ≤ B ≤ C, the longest axis of this triaxial ellipsoid coincides with the ĉ direction (all the angular momentum is along the axis of largest principal moment of inertia) and the shortest axis of the energy ellipsoid coincides with the â axis (all the angular momentum is along the smallest moment of inertia). Without actually solving the Lagrange equations, we have found strong constraints on the evolution of the components of the angular momentum on the principal axes.

To determine how the system evolves along these intersection curves we have to use the equations of motion. We observe that the evolution of the components of the angular momentum on the principal axes depends only on the components of the angular momentum on the principal axes, even though the values of these components are not enough to completely specify the dynamical state. Apparently the dynamics of these components is self-contained, and we will see that it can be described in terms of a set of differential equations whose only dynamical variables are the components of the angular momentum on the principal axes (see section 2.9).

We note that there are two axes for which the intersection curves shrink to a point if we hold the energy constant and vary the magnitude of the angular momentum. If the angular momentum starts at these points, the conserved quantities constrain the angular momentum to stay there. These points are equilibrium points for the body components of the angular momentum. However, they are not equilibrium points for the system as a whole. At these points the body is still rotating even though the body components of the angular momentum are not changing. This kind of equilibrium is called a relative equilibrium. We can also see that if the angular momentum is initially slightly displaced from one of these relative equilibria, then the angular momentum is constrained to stay near it on one of the intersection curves. The angular momentum vector is fixed in space, so the principal axis of the equilibrium point of the body rotates stably about the angular momentum vector.

At the principal axis with intermediate moment of inertia, the $\hat{b}$ axis, the intersection curves appear to cross. As we observed, the dynamics of the components of the angular momentum on the principal axes forms a self-contained dynamical system. Trajectories of a dynamical system cannot cross,⁹ so the most that can happen is that if the equations of motion carry the system along the intersection curve then the system can approach the crossing point only asymptotically. So without solving any equations we can deduce that the point of crossing is another relative equilibrium. If the angular momentum is initially aligned with the intermediate axis, then it stays aligned. If the system is slightly displaced from the intermediate axis, then the evolution along the intersection curve will take the system far from the relative equilibrium. So rotation about the axis of intermediate moment of inertia is unstable—initial displacements of the angular momentum, however small initially, become large. Again, the angular momentum vector is fixed in space, but now the principal axis with the intermediate principal moment does not stay close to the angular momentum, so the body executes a complicated tumbling motion.

This gives some insight into the mystery of the thrown book mentioned at the beginning of the chapter. If one throws a book so that it is initially rotating about either the axis with the largest moment of inertia or the axis with the smallest moment of inertia (the shortest and longest physical axes, respectively), the book rotates regularly about that axis. However, if the book is thrown so that it is initially rotating about the axis of intermediate moment of inertia (the intermediate physical axis), then it tumbles, however carefully it is thrown. You can try it with this book (but put a rubber band or string around it first).

Before moving on, we can make some further physical deductions. Suppose a freely rotating body is subject to some sort of internal friction that dissipates energy but conserves the angular momentum. For example, real bodies flex as they spin. If the spin axis moves with respect to the body then the flexing changes with time, and this changing distortion converts kinetic energy of rotation into heat. Internal processes do not change the total angular momentum of the system. If we hold the magnitude of the angular momentum fixed but gradually decrease the energy, then the curve of intersection on which the system moves gradually deforms. For a given angular momentum there is a lower limit on the energy: the energy cannot be so low that there are no intersections. For this lowest energy the intersection of the angular momentum sphere and the energy ellipsoid is a pair of points on the axis of maximum moment of inertia. With energy dissipation, a freely rotating physical body eventually ends up with the lowest energy consistent with the given angular momentum, which is rotation about the principal axis with the largest moment of inertia (typically the shortest physical axis).

Thus, we expect that given enough time all freely rotating physical bodies will end up rotating about the axis of largest moment of inertia. You can demonstrate this to your satisfaction by twirling a small bottle containing some viscous fluid, such as correction fluid. What you will find is that, whatever spin you try to put on the bottle, it will reorient itself so that the axis of the largest moment of inertia is aligned with the spin axis. Remarkably, this is very nearly true of almost every body in the solar system for which there is enough information to decide. The deviations from principal axis rotation for the Earth are tiny: the angle between the angular momentum vector and the ĉ axis for the Earth is less than one arc-second.¹⁰ In fact, the evidence is that all of the planets, the Moon and all the other natural satellites, and almost all of the asteroids rotate very nearly about the largest moment of inertia. We have deduced that this is to be expected using an elementary argument. There are exceptions. Comets typically do not rotate about the largest moment. As they are heated by the sun, material spews out from localized jets, and the back reaction from these jets changes the rotation state. Among the natural satellites, the only known exception is Saturn's satellite Hyperion, which is tumbling chaotically. Hyperion is especially out of round and subject to strong gravitational torques from Saturn.

2.9   Euler's Equations

For a free rigid body we have seen that the components of the angular momentum on the principal axes comprise a self-contained dynamical system: the variation of the principal axis components depends only on the principal axis components. Here we derive equations that govern the evolution of these components.

The starting point for the derivation is the conservation of the vector angular momentum. The components of the angular momentum on the principal axes are¹¹

$$\begin{array}{cc}\overline{L}\prime =I\prime \mathit{\omega }\prime ,& \left(2.61\right)\end{array}$$

where ω′ is composed of the components of the angular velocity vector on the principal axes and I′ is the matrix representation of the inertia tensor with respect to the principal axis basis:

$$\begin{array}{cc}I\prime =\left(\begin{array}{ccc}A& 0& 0\\ 0& B& 0\\ 0& 0& C\end{array}\right).& \left(2.62\right)\end{array}$$

The body components of the angular momentum L′ are related to the components L on the fixed rectangular basis ê_i by

$$\begin{array}{cc}\overline{L}=\mathbf{M}\overline{L}\prime ,& \left(2.63\right)\end{array}$$

where M is the matrix representation of the rotation that carries the body and all vectors attached to the body from the reference orientation of the body to the actual orientation.

The vector angular momentum is conserved for free rigid-body motion, and so are its components on a fixed rectangular basis. So, along solution paths,

$$\begin{array}{cc}0=D\overline{L}=DM\overline{L}\prime +MD\overline{L}\prime .& \left(2.64\right)\end{array}$$

Solving, we find

$$\begin{array}{cc}D\overline{L}\prime =-M^{\mathcal{T}}DM\overline{L}\prime .& \left(2.65\right)\end{array}$$

In terms of ω′ this is

$$\begin{array}{lll}I\prime D\mathit{\omega }\prime \hfill & =-M^{\mathcal{T}}DMI\prime \mathit{\omega }\prime \hfill & \hfill \\ \hfill & =-M^{\mathcal{T}}\mathcal{A}\left(M\mathit{\omega }\prime \right)MI\prime \mathit{\omega }\prime ,\hfill & \left(2.66\right)\hfill \end{array}$$

where we have used equation (2.21) to write DM in terms of $\mathcal{A}$. The function $\mathcal{A}$ has the property¹²

$$\begin{array}{cc}{R}^{\mathcal{T}}\mathcal{A}\left(Rv\right)R=\mathcal{A}\left(v\right)& \left(2.67\right)\end{array}$$

for any vector with components v and any rotation with matrix representation R. Using this property of $\mathcal{A}$, we find Euler's equations:

$$\begin{array}{cc}I\prime D\mathit{\omega }\prime =-\mathcal{A}(\mathit{\omega }\prime )I\prime \mathit{\omega }\prime .& \left(2.68\right)\end{array}$$

Euler's equations give the time derivative of the body components of the angular velocity vector entirely in terms of the angular velocity components and the principal moments of inertia. Let ω^a, ω^b, and ω^c denote the components of the angular velocity vector on the principal axes. Then Euler's equations can be written as the component equations

$$\begin{array}{ll}AD{\omega }^a=\left(B-C\right){\omega }^b{w}^c\hfill & \hfill \\ BD{\omega }^b=\left(C-A\right){\omega }^c{\omega }^a\hfill & \hfill \\ CD{\omega }^c=\left(A-B\right){\omega }^a{\omega }^b.\hfill & \left(2.69\right)\hfill \end{array}$$

Alternatively, we can rewrite Euler's equations in terms of the components of the angular momentum on the principal axes

$$\begin{array}{ll}D{L}_a=\left(\frac{1}{C}-\frac{1}{B}\right)L_b{L}_c\hfill & \hfill \\ D{L}_b=\left(\frac{1}{A}-\frac{1}{C}\right)L_a{L}_c\hfill & \hfill \\ D{L}_a=\left(\frac{1}{B}-\frac{1}{A}\right)L_a{L}_b.\hfill & \left(2.70\right)\hfill \end{array}$$

These equations confirm that the time derivatives of the components of the angular momentum on the principal axes depend only on the components of the angular momentum on the principal axes.

Euler's equations are very simple, but they do not completely determine the evolution of a rigid body—they do not give the spatial orientation of the body. However, equation (2.21) and property (2.67) can be used to relate the derivative of the orientation matrix to the body components of the angular velocity vector:

$$\begin{array}{cc}DM=M\mathcal{A}(\mathit{\omega }\prime ).& \left(2.71\right)\end{array}$$

A straightforward method of using these equations is to integrate them componentwise as a set of nine first-order ordinary differential equations, with initial conditions determining the initial configuration matrix. Together with Euler's equations, which describe how the body components of the angular velocity vector change with time, this system of equations governing the motion of a rigid body is complete. However, the reader will no doubt have noticed that this approach is rather wasteful. The fact that the orientation matrix can be specified with only three parameters has not been taken into account. We should be integrating three equations for the orientation, given ω′, not nine. To accomplish this we once again need to parameterize the configuration matrix.

For example, we can use Euler angles to parameterize the orientation:

$$\begin{array}{cc}\mathcal{M}\left(\theta ,\phi ,\psi \right)=R_z\left(\phi \right)R_x\left(\theta \right)R_z\left(\psi \right).& \left(2.72\right)\end{array}$$

We form M by composing $\mathcal{M}$ with an Euler coordinate path. Equation (2.71) can then be used to solve for Dθ, Dφ, and Dψ. We find

$$\begin{array}{cc}\left(\begin{array}{l}D\theta \hfill \\ D\phi \hfill \\ D\psi \hfill \end{array}\right)=\frac{1}{\sin \theta }\left(\begin{array}{ccc}\cos \psi \sin \theta & -\sin \psi \sin \theta & 0\\ \sin \psi & \cos \psi & 0\\ -\sin \psi \cos \theta & -\cos \psi \cos \theta & \sin \theta \end{array}\right)\left(\begin{array}{c}{\omega }^a\\ {\omega }^b\\ {\omega }^c\end{array}\right).& \left(2.73\right)\end{array}$$

This gives us the desired equation for the orientation. Note that it is singular for θ = 0, as are Lagrange's equations. So Euler's equations using Euler angles for the configuration have the same problem as did the Lagrange equations using Euler angles. Again, this is a manifestation of the fact that for θ = 0 the orientation depends only on φ+ψ. The singularity in the equations of motion for θ = 0 does not correspond to anything funny in the motion of the rigid body. A practical solution to the singularity problem is to choose another set of Euler-like angles that have a singularity in a different place, and switch from one to the other when the going gets tough.

Exercise 2.12:

Fill in the details of the derivation of equation (2.73). You may want to use the computer to help with the algebra.

Euler's equations for forced rigid bodies

Euler's equations were derived for a free rigid body. In general, we must be able to deal with external forcing. How do we do this? First, we derive expressions for the vector torque. Then we include the vector torque in the Euler equations.

We derive the vector torque in a manner analogous to the derivation of the vector angular momentum. That is, we derive one component and then argue that since the coordinate system is arbitrary, all components have the same form.

Suppose we have a rigid body subject to some potential energy that depends only on time and the configuration. A Lagrangian is L = T − V. If we use the Euler angles as generalized coordinates, the last of the three active Euler rotations that define the orientation is a rotation about the ẑ axis by the angle φ. The Lagrange equation for φ gives¹³

$$\begin{array}{cc}D{p}_{\phi }(t)=-{\partial }_{1,1}V(t;\theta (t),\phi (t),\psi (t)).& (2.74)\end{array}$$

If we define T_z, the component of the torque about the z axis, to be minus the derivative of the potential energy with respect to the angle of rotation of the body about the z axis,

$$\begin{array}{cc}{T}_z\left(t\right)=-{\partial }_{1,1}V\left(t;\theta \left(t\right),\phi \left(t\right),\psi \left(t\right)\right),& \left(2.75\right)\end{array}$$

then we see that

$$\begin{array}{cc}D{p}_{\phi }\left(t\right)=T_z\left(t\right).& \left(2.76\right)\end{array}$$

We have already identified the momentum conjugate to φ as one component, L_z, of the vector angular momentum $\overrightarrow{L}$ (see section 2.8), so

$$\begin{array}{cc}D{L}_z\left(t\right)=T_z.& \left(2.77\right)\end{array}$$

Since the orientation of the reference rectangular basis vectors is arbitrary, we may choose them any way that we please. Thus if we want any component of the vector torque, we may choose the z-axis so that we can compute it in this way. We can conclude that the vector torque gives the rate of change of the vector angular momentum

$$\begin{array}{cc}D\overrightarrow{L}=\overrightarrow{T}.& \left(2.78\right)\end{array}$$

Having obtained a general prescription for the vector torque, we address how the vector torque may be included in Euler's equations. Euler's equations expressed the fact that the vector angular momentum is conserved. Let's return to that calculation, but now include a torque with components T arranged as a column matrix:

$$\begin{array}{cc}D\overline{L}=\overline{T}=DM\overline{L}\prime +MD\overline{L}\prime .& \left(2.79\right)\end{array}$$

Carrying out the same steps as before, we find

$$\begin{array}{ll}{T}_a=D{L}_a-\left(\frac{1}{C}-\frac{1}{B}\right)L_b{L}_c\hfill & \hfill \\ T_b=D{L}_b-\left(\frac{1}{A}-\frac{1}{C}\right)L_a{L}_c\hfill & \hfill \\ T_c=D{L}_a-\left(\frac{1}{B}-\frac{1}{A}\right)L_a{L}_b,\hfill & \left(2.80\right)\hfill \end{array}$$

where the components of the torque on the principal axes are

$$\begin{array}{cc}\overline{T}\prime =M^{-1}\overline{T}.& \left(2.81\right)\end{array}$$

In terms of ω′ this is

$$\begin{array}{cc}I\prime D\mathit{\omega }\prime +\mathcal{A}(\mathit{\omega }\prime )I\prime \mathit{\omega }\prime =\overline{T}\prime ;& \left(2.82\right)\end{array}$$

in components,

$$\begin{array}{ll}AD{\omega }^a-\left(B-C\right){\omega }^b{\omega }^c=T_a\hfill & \hfill \\ BD{\omega }^b-\left(C-A\right){\omega }^c{\omega }^a=T_b\hfill & \hfill \\ CD{\omega }^c-\left(A-B\right){\omega }^a{\omega }^b=T_c.\hfill & \left(2.83\right)\hfill \end{array}$$

Note that the torque entered only the equations for the body angular momentum and for the body angular velocity vector. The equations that relate the derivative of the orientation to the angular velocity vector are not modified by the torque. In a sense, Euler's equations contain the dynamics, and the equations governing the orientation are kinematic. Of course, Lagrange's equations must be modified by the potential that gives rise to the torques; in this sense Lagrange's equations contain both dynamics and kinematics.

Exercise 2.13: Bicycle wheel

a. Imagine that you are holding a bicycle wheel by the axle (in both hands) and the wheel is spinning so that the top edge is going away from your face. If you torque the wheel by pushing down with your right hand and pulling up with your left hand the wheel will precess. Which way does it try to turn?

b. A free bicycle wheel rolls on a horizontal surface. If it starts to tilt, the torque from gravity will cause the wheel to turn. Which way will it turn? The reasoning that applied to part a does not directly apply to the rolling bicycle wheel, which is not a holonomic system. However, it is interesting to think about whether the behavior of the two systems is related.

2.10   Axisymmetric Tops

We have all played with a top at one time or another. For the purposes of analysis we will consider an idealized top that does not wander around. Thus, an ideal top is a rotating rigid body, one point of which is fixed in space. Furthermore, the center of mass of the top is not at the fixed point, which is the center of rotation, and there is a uniform gravitational acceleration.

For our top we can take the Lagrangian to be the difference of the kinetic energy and the potential energy. We already know how to write the kinetic energy—what is new here is that we must express the potential energy in terms of the configuration. In the case of a body in a uniform gravitational field this is easy. The potential energy is the sum of “mgh” for all the constituent particles:

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{m}_{\alpha }g{h}_{\alpha },}& \left(2.84\right)\end{array}$$

where g is the gravitational acceleration, h_α = ${\overrightarrow{x}}_{\alpha }$ · ẑ, and the unit vector ẑ indicates which way is up. Rewriting the vector to the constituents in terms of the vector $\overrightarrow{X}$ to the center of mass, the potential energy is

$$\begin{array}{ll}{\displaystyle \sum _{\alpha }{m}_{\alpha }g\left(\overrightarrow{X}+{\overrightarrow{\xi }}_{\alpha }\right)·\hat{z}}\hfill & \hfill \\ =gM\overrightarrow{X}·\hat{z}+g\left({\displaystyle \sum _{\alpha }{m}_{\alpha }{\overrightarrow{\xi }}_{\alpha }}\right)·\hat{z}\hfill & \hfill \\ =gM\overrightarrow{X}·\hat{z},\hfill & \left(2.85\right)\hfill \end{array}$$

where the last sum is zero because the center of mass is the origin of ${\overrightarrow{\xi }}_{\alpha }$. So the potential energy of a body in a gravitational field with uniform acceleration is very simple: it is just M gh, where M is the total mass and $h=\overrightarrow{X}·\hat{z}$ is the height of the center of mass.

Here we consider an axisymmetric top (see figure 2.4). Such a top has an axis of symmetry of the mass distribution, so the center of mass is on the symmetry axis and the fixed point is also on the axis of symmetry.

In order to write the Lagrangian we need to choose a set of generalized coordinates. If we choose them well we can take advantage of the symmetries of the problem. If the Lagrangian does not depend on a particular coordinate, the conjugate momentum is conserved, and the complexity of the system is reduced.

The axisymmetric top has two apparent symmetries. The fact that the mass distribution is axisymmetric implies that neither the kinetic nor the potential energy is sensitive to the orientation of the top about that symmetry axis. Additionally, the kinetic and potential energy are insensitive to a rotation of the physical system about the vertical axis, because the gravitational field is uniform.

We can take advantage of these symmetries by choosing appropriate coordinates, and we already have a coordinate system that does the job—the Euler angles.¹⁴ We choose the reference orientation so that the symmetry axis is vertical. The first Euler angle, ψ, expresses a rotation about the symmetry axis. The next Euler angle, θ, is the tilt of the symmetry axis of the top from the vertical. The third Euler angle, φ, expresses a rotation of the top about the ẑ axis. The symmetries of the problem imply that the first and third Euler angles do not appear in the Lagrangian. As a consequence the momenta conjugate to these angles are conserved quantities. Let's work out the details.

First, we develop the Lagrangian explicitly. The general form of the kinetic energy about a fixed point is given by equation 2.41. The top is constrained so that it pivots about a fixed point that is not at the center of mass. So the moments of inertia that enter the kinetic energy are the moments of inertia of the top with respect to the pivot point, not with respect to the center of mass. If we know the moments of inertia about the center of mass we can write the moments of inertia about the pivot in terms of them (see exercise 2.2 on Steiner's theorem). So let's assume the principal moments of inertia of the top about the pivot are A, B, and C, and A = B because of the symmetry.¹⁵ We can use the computer to help us figure out the Lagrangian for this special case:

(show-expression
  ((T-body-Euler 'A 'A 'C)
   (up 't
       (up 'theta 'phi 'psi)
       (up 'thetadot 'phidot 'psidot))))

We can rearrange this a bit to get

$$\begin{array}{ll}T\left(t;\theta ,\phi ,\psi ;\dot{\theta },\dot{\phi },\dot{\psi }\right)\hfill & \hfill \\ =\frac{1}{2}A\left({\dot{\theta }}^2+{\dot{\phi }}^2{\sin }^2\theta \right)+\frac{1}{2}C{\left(\dot{\psi }+\dot{\phi }\cos \theta \right)}^2.\hfill & \left(2.86\right)\hfill \end{array}$$

In terms of Euler angles, the potential energy is

$$\begin{array}{cc}V\left(t;\theta ,\phi ,\psi ;\dot{\theta },\dot{\phi },\dot{\psi }\right)=MgR\cos \theta ,& \left(2.87\right)\end{array}$$

where R is the distance of the center of mass from the pivot. The Lagrangian is L = T − V. We see that the Lagrangian is indeed independent of ψ and φ, as expected.

There is no particular reason to look at the Lagrange equations. We can assign that job to the computer when needed. However, we have already seen that it may be useful to examine the conserved quantities associated with the symmetries.

The energy is conserved, because the Lagrangian has no explicit time dependence. Also, the energy is the sum of the kinetic and potential energy E = T + V, because the kinetic energy is a homogeneous quadratic form in the generalized velocities. The energy is

$$\begin{array}{cc}E=\frac{1}{2}A\left({\dot{\theta }}^2+{\dot{\phi }}^2{\sin }^2\theta \right)+\frac{1}{2}C{\left(\dot{\psi }+\dot{\phi }\cos \theta \right)}^2+MgR\cos \theta .& \left(2.88\right)\end{array}$$

Two of the generalized coordinates do not appear in the Lagrangian, so there are two conserved momenta. The momentum conjugate to φ is

$$\begin{array}{cc}{p}_{\phi }=\left(A{\left(\sin \theta \right)}^2+C{\left(\cos \theta \right)}^2\right)\dot{\phi }+C\dot{\psi }\cos \theta .& \left(2.89\right)\end{array}$$

The momentum conjugate to ψ is

$$\begin{array}{cc}{p}_{\psi }=C\left(\dot{\psi }+\dot{\phi }\cos \theta \right).& \left(2.90\right)\end{array}$$