Mass and Inertia#
Learning Objectives#
After reading this chapter the reader should know:
What the inertia matrix is and how it is constructed from the mass distribution of a body.
How the inertia matrix transforms when changing reference frames or reference points.
What angular momentum is and how it relates to the inertia matrix and angular velocity.
The Newton-Euler equations of motion for a rigid body, both about the center of mass and about an arbitrary point.
How to apply the parallel axis theorem.
Introduction#
Picture a figure skater doing a spin. As they pull their arms in close to their body, they spin noticeably faster. The distribution of mass is controlling how easy or hard it is to change the rotation, and that is represented by the inertia matrix. Just as mass is the scalar quantity that resists changes to linear velocity (\(F = ma\)), the inertia matrix is the 3×3 tensor quantity that resists changes to angular velocity. When the skater pulls in their arms, the mass moves closer to the axis of rotation and the inertia decreases, so the angular velocity must increase to conserve angular momentum.
Given a rotating body, what inertia does it present to an applied torque? The answer to this question is what drives the rotational equations of motion (Newton-Euler), which we also derive in this chapter. By the end of the page you will be able to write down and compute the full equations of motion for a rigid body in 3D.
A spinning top#
Throughout this chapter we use a spinning top (often just referred to as “top”), a solid cone with a short cylindrical stem, as our running example, building up its inertia properties step by step.
import numpy as np
# Cone (body)
m_cone = 0.4 # kg
R_cone = 0.04 # m
H_cone = 0.06 # m
# Stem (handle)
m_stem = 0.1 # kg
R_stem = 0.008 # m
L_stem = 0.04 # m
# Bolt (used in some examples to break symmetry)
m_bolt = 0.02 # kg
r_bolt = np.array([R_cone, 0, H_cone]) # m
Mass and Center of Mass#
Before tackling the full inertia matrix, let us recall the simpler scalar quantities. The total mass of a body is obtained by integrating the mass density over the body’s volume:
The center of mass is the mass-weighted average position:
For a system of \(N\) discrete point masses \(m_k\) at positions \(\mathbf{r}_k\), these integrals become sums:
Our spinning top is a solid cone (the body) with a short cylindrical stem (the handle). The cone’s tip rests on the surface at the origin, and the spin axis points along \(z\). Both components are axisymmetric, so their centers of mass lie on the \(z\)-axis. We can find the combined center of mass by treating each component as a point mass located at its own center of mass.
# Centers of mass of each component
z_cm_cone = 3 * H_cone / 4
r_cm_cone = np.array([0, 0, z_cm_cone])
z_cm_stem = H_cone + L_stem / 2
r_cm_stem = np.array([0, 0, z_cm_stem])
masses = np.array([m_cone, m_stem])
positions = np.array([r_cm_cone, r_cm_stem])
m_total = np.sum(masses)
r_cm_top = np.average(positions, axis=0, weights=masses)
z_cm_top = r_cm_top[2]
print(f"Total mass: {m_total:.2f} kg")
print(f"Cone CoM (above tip): z = {z_cm_cone:.4f} m")
print(f"Stem CoM (above tip): z = {z_cm_stem:.4f} m")
print(f"Combined CoM (above tip): z = {z_cm_top:.4f} m")
Total mass: 0.50 kg
Cone CoM (above tip): z = 0.0450 m
Stem CoM (above tip): z = 0.0800 m
Combined CoM (above tip): z = 0.0520 m
Exercise
A small mounting bolt of mass is attached at the rim of the cone at its top.
print(f"Bolt mass (kg): {m_bolt}")
print(f"Bolt position (m): {r_bolt}")
Bolt mass (kg): 0.02
Bolt position (m): [0.04 0. 0.06]
Find the center of mass of the combined system (cone + stem + bolt). By how many millimetres does the bolt shift the CoM in the \(x\)-direction?
Solution
masses = np.array([m_cone, m_stem, m_bolt])
positions = np.array([r_cm_cone, r_cm_stem, r_bolt])
r_cm_boltedtop = np.average(positions, axis=0, weights=masses)
print(f"Center of mass: {r_cm_boltedtop}")
print(f"x-shift due to bolt: {r_cm_boltedtop[0]*1000:.3f} mm")
Center of mass: [0.00153846 0. 0.05230769]
x-shift due to bolt: 1.538 mm
The Inertia Matrix#
Definition via Integral#
For rotation, the role of scalar mass is played by the inertia matrix \(\mathbf{M}_{b/c}\), which encodes how the mass of body \(b\) is distributed around its center of mass \(c\).
Inertia Matrix
The inertia matrix of body \(b\) about center of mass \(c\), expressed in frame \(i\), is
where \(\mathbf{r}^i\) is the position of the mass element relative to \(c\) expressed in frame \(i\), and \([\mathbf{r}]^\times\) is the skew-symmetric cross-product matrix such that \([\mathbf{r}]^\times \mathbf{v} = \mathbf{r} \times \mathbf{v}\) for any vector \(\mathbf{v}\).
The skew-symmetric matrix of a vector \(\mathbf{r} = [x, y, z]^\top\) is:
Let us build this in SymPy and compute the inertia contribution of a single point mass at a general position \((x, y, z)\):
import sympy as sm
sm.init_printing(use_latex='mathjax')
x, y, z, m = sm.symbols('x y z m', real=True)
r = sm.Matrix([x, y, z])
M_point = m * (r.dot(r) * sm.eye(3) - r * r.T)
M_point
We can verify this equals the skew-matrix formula:
def skew_sm(v):
return sm.Matrix([
[ 0, -v[2], v[1]],
[ v[2], 0, -v[0]],
[-v[1], v[0], 0],
])
def skew_np(v):
return np.array([
[ 0, -v[2], v[1]],
[ v[2], 0, -v[0]],
[-v[1], v[0], 0]
])
r_cross = skew_sm(r)
M_from_skew = -m * r_cross * r_cross
sm.simplify(M_from_skew)
General Properties#
The inertia matrix is always symmetric (follows directly from the definition) and positive definite, meaning the rotational kinetic energy is always strictly positive:
Positive definiteness guarantees all eigenvalues of \(\mathbf{M}\) are strictly positive. Let’s ensure that this is the case for our spinning top:
# Cone and stem inertia formulas
Izz_cone = 3 * m_cone * R_cone**2 / 10
Ixx_cone = Iyy_cone = 3 * m_cone * (4*R_cone**2 + H_cone**2) / 80
Izz_stem = m_stem * R_stem**2 / 2
Ixx_stem = Iyy_stem = m_stem * (3*R_stem**2 + L_stem**2) / 12
# Body-frame inertia for the cone and stem
M_cone = np.diag([Ixx_cone, Iyy_cone, Izz_cone])
M_stem = np.diag([Ixx_stem, Iyy_stem, Izz_stem])
# Vectors from the top's cm to each component's cm
r_top_to_cone = r_cm_cone - r_cm_top
r_top_to_stem = r_cm_stem - r_cm_top
# Parallel axis theorem to calculate the total top inertia matrix
# The theorem will be discussed in greater detail later
M_top = (
M_cone - m_cone * skew_np(r_top_to_cone) @ skew_np(r_top_to_cone) +
M_stem - m_stem * skew_np(r_top_to_stem) @ skew_np(r_top_to_stem)
)
print("Top inertia in the body frame (kg·m²):")
print(np.round(M_top * 1e6, 2), " (×10⁻⁶)")
Top inertia in the body frame (kg·m²):
[[262.93 0. 0. ]
[ 0. 262.93 0. ]
[ 0. 0. 195.2 ]] (×10⁻⁶)
Frame Transformations#
The inertia matrix depends on which frame the position vectors \(\mathbf{r}\) are expressed in. The body frame \(b\) is usually the most convenient for computation (the geometry is fixed), but the Newton-Euler equations are defined in the inertial frame \(i\). The transformation between the two is:
where \(\mathbf{R}^i_b\) is the rotation matrix from frame \(b\) to frame \(i\), and \(\mathbf{R}^b_i = (\mathbf{R}^i_b)^\top\).
A rotation never changes the physical mass distribution, so the eigenvalues of the inertia matrix are frame-independent. Let us verify this by comparing the inertia of the spinning top’s body in the body frame with its inertia in the inertial frame after being tilted 30° about the y-axis:
# Tilt the top 30° about the y-axis
theta = np.radians(30)
Ry = np.array([
[ np.cos(theta), 0.0, np.sin(theta)],
[ 0.0, 1.0, 0.0],
[-np.sin(theta), 0.0, np.cos(theta)],
])
M_top_inertial = Ry @ M_top @ Ry.T
print("\nInertia in inertial frame after 30° tilt (kg·m²):")
print(np.round(M_top_inertial * 1e6, 2), " (×10⁻⁶)")
eigs_body = np.sort(np.linalg.eigvalsh(M_top))
eigs_inertial = np.sort(np.linalg.eigvalsh(M_top_inertial))
print(f"\nEigenvalues body frame: {np.round(eigs_body*1e6, 4)} (×10⁻⁶)")
print(f"Eigenvalues inertial frame: {np.round(eigs_inertial*1e6, 4)} (×10⁻⁶)")
Inertia in inertial frame after 30° tilt (kg·m²):
[[246. 0. -29.33]
[ 0. 262.93 0. ]
[-29.33 0. 212.13]] (×10⁻⁶)
Eigenvalues body frame: [195.2 262.9333 262.9333] (×10⁻⁶)
Eigenvalues inertial frame: [195.2 262.9333 262.9333] (×10⁻⁶)
The Parallel Axis Theorem#
You computed the inertia matrix about the center of mass, but you need it about the wheel hub, a joint, or the end of a rod. The parallel axis theorem lets you shift the reference point without re-integrating.
Parallel Axis Theorem
Given the inertia matrix \(\mathbf{M}_{b/c}\) about the center of mass \(c\), the inertia matrix about a new point \(o\) is
where \(\mathbf{r}_{o/c}\) is the vector from \(c\) to \(o\).
In 2D, rotating about a fixed axis, the theorem reduces to the scalar form:
where \(d\) is the perpendicular distance between the original axis (through \(c\)) and the new axis (through \(o\)). The inertia always increases when moving away from the center of mass, which means the correction \(-m[\mathbf{r}]^\times[\mathbf{r}]^\times\) is positive semi-definite.
Angular Momentum#
About the Center of Mass#
Angular momentum is the rotational analogue of linear momentum \(\mathbf{p} = m\mathbf{v}\). The angular momentum of body \(b\) about its center of mass \(c\) is:
Note that \(\mathbf{h}_{b/c}\) is generally not parallel to \(\boldsymbol{\omega}_{b/i}\) unless the rotation happens to be about a principal axis.
For our axisymmetric spinning top (no bolt), the inertia matrix in the body frame is diagonal. Spinning along the \(z\)-axis therefore gives angular momentum exactly parallel to \(\boldsymbol{\omega}\):
omega = np.array([0.0, 0.0, 100.0])
h = M_top @ omega
cos_angle = np.dot(h, omega) / (np.linalg.norm(h) * np.linalg.norm(omega))
angle_deg = np.degrees(np.arccos(np.clip(cos_angle, -1.0, 1.0)))
print("--- Balanced cone (no bolt) ---")
print(f"omega: {omega} rad/s")
print(f"h: {np.round(h, 6)} kg·m²/s")
print(f"Angle between h and omega: {angle_deg:.4f}°")
--- Balanced cone (no bolt) ---
omega: [ 0. 0. 100.] rad/s
h: [0. 0. 0.01952] kg·m²/s
Angle between h and omega: 0.0000°
Now attach the bolt at the rim. The bolt’s inertia contribution has off-diagonal terms because its position \([R\_cone, 0, H\_cone]\) is not on the spin axis:
# Vectors from the bolted top's cm to each component's cm
r_boltedtop_to_cone = r_cm_cone - r_cm_boltedtop
r_boltedtop_to_stem = r_cm_stem - r_cm_boltedtop
r_boltedtop_to_bolt = r_bolt - r_cm_boltedtop
# Parallel axis theorem
M_boltedtop = (
M_cone - m_cone * skew_np(r_boltedtop_to_cone) @ skew_np(r_boltedtop_to_cone) +
M_stem - m_stem * skew_np(r_boltedtop_to_stem) @ skew_np(r_boltedtop_to_stem) -
m_bolt * skew_np(r_boltedtop_to_bolt) @ skew_np(r_boltedtop_to_bolt)
)
h_bolted = M_boltedtop @ omega
cos_angle = np.dot(h_bolted, omega) / (np.linalg.norm(h_bolted) * np.linalg.norm(omega))
angle_deg = np.degrees(np.arccos(np.clip(cos_angle, -1.0, 1.0)))
print("--- Unbalanced cone (with bolt) ---")
print(f"omega: {omega} rad/s")
print(f"h: {np.round(h_bolted, 6)} kg·m²/s")
print(f"Angle between h and omega: {angle_deg:.4f}°")
--- Unbalanced cone (with bolt) ---
omega: [ 0. 0. 100.] rad/s
h: [-0.000615 0. 0.022597] kg·m²/s
Angle between h and omega: 1.5600°
The bolt shifts the angular momentum away from the spin axis. This causes a wobble when spinning the top with a bolt on its rim.
About a Different Point#
When computing angular momentum about a point \(o\) that is not the center of mass, an extra term appears due to the translational motion of the center of mass:
where \(\mathbf{r}_{c/o}\) is the position of \(c\) relative to \(o\), and \(\mathbf{v}_{c/i}\) is the velocity of the center of mass in the inertial frame.
Time Derivative#
Taking the time derivative of angular momentum and setting it equal to the applied torque gives:
The second term, \(\boldsymbol{\omega} \times (\mathbf{M}\boldsymbol{\omega})\), vanishes when \(\boldsymbol{\omega}\) is parallel to \(\mathbf{M}\boldsymbol{\omega}\). This is the case when spinning about a principal axis.
Newton-Euler Equations of Motion#
Now we can write down the full equations of motion for a rigid body. There are two vector equations. One for translation and one for rotation.
About the Center of Mass#
Newton-Euler Equations (About Center of Mass)
The translational equation is simply \(F = ma\) applied to the center of mass \(c\). In block matrix form both equations read:
About a Different Point o#
When the body is constrained to rotate about a fixed point \(o\) that is not the center of mass, the equations change because the offset \(\mathbf{r}^i_{c/o}\) between the center of mass and the pivot couples translation and rotation:
In block matrix form:
Simulation: Spinning Top#
The equations of motion are most instructive in motion. With the tip fixed at the origin, \(\mathbf{a}_o = \mathbf{0}\), and the torque equation from (58) reduces to
where the only applied torque is gravity acting at the center of mass: ..
To integrate this forward in time we track the orientation as a unit quaternion \(\mathbf{q} = [q_w,\, q_x,\, q_y,\, q_z]^\top\) whose kinematics are
The ODE state is \(\mathbf{y} = [q_w,\, q_x,\, q_y,\, q_z,\, \omega_x,\, \omega_y,\, \omega_z]^\top\). We start from zero rotation (perfectly upright) and a 200 rad/s spin.
from scipy.integrate import solve_ivp
from scipy.spatial.transform import Rotation
g = 9.81
g_world = np.array([0.0, 0.0, -g])
# Inertia about tip using the parallel axis theorem
M_tip = (
M_cone - m_cone * skew_np(r_cm_cone) @ skew_np(r_cm_cone) +
M_stem - m_stem * skew_np(r_cm_stem) @ skew_np(r_cm_stem)
)
r_cm_body = np.array([0.0, 0.0, z_cm_top])
def quat_mult(p, r):
"""Quaternion product of p and r, both in [w, x, y, z] format."""
pw, px, py, pz = p
rw, rx, ry, rz = r
return np.array([
pw*rw - px*rx - py*ry - pz*rz,
pw*rx + px*rw + py*rz - pz*ry,
pw*ry - px*rz + py*rw + pz*rx,
pw*rz + px*ry - py*rx + pz*rw,
])
def top_ode(t, y, M, M_inv, r_cm, m):
q = y[:4]
omega = y[4:]
R = Rotation.from_quat([*q[1:], q[0]]).as_matrix()
tau_body = R.T @ np.cross(R @ r_cm, m * g_world)
alpha = M_inv @ (tau_body - np.cross(omega, M @ omega))
dq = 0.5 * quat_mult(q, np.array([0.0, *omega]))
return [*dq, *alpha]
# Initial conditions: zero rotation, 200 rad/s spin about body-z
y0 = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 200.0]
t_span = (0.0, 4.0)
t_eval = np.linspace(0.0, 4.0, 480)
sol_top = solve_ivp(
top_ode, t_span, y0, t_eval=t_eval, method='RK45',
args=(M_tip, np.linalg.inv(M_tip), r_cm_body, m_total),
rtol=1e-8, atol=1e-10,
)
quats_top = sol_to_yup_quats(sol_top)
renderer_top, action_top = make_top_animation(quats_top, sol_top.t)
renderer_top
With a perfectly balanced top starting from zero rotation, the spin axis is aligned with gravity and with the principal axis, and the top spins without any wobble.
Now attach the bolt. The initial rotation is still zero (perfectly upright, 200 rad/s spin), but the bolt shifts the center of mass off the symmetry axis. Gravity now acts at a point that is no longer on the spin axis, immediately creating a torque. The result is visible wobble, even though the starting orientation is identical to the balanced case. This is the signature of a broken principal axis:
# Bolt's point-mass inertia contribution about the tip
M_bolt_tip = -m_bolt * skew_np(r_bolt) @ skew_np(r_bolt)
M_tip_bolted = M_tip + M_bolt_tip
m_bolted = m_total + m_bolt
sol_bolted = solve_ivp(
top_ode, t_span, y0, t_eval=t_eval, method='RK45',
args=(M_tip_bolted, np.linalg.inv(M_tip_bolted), r_cm_boltedtop, m_bolted),
rtol=1e-8, atol=1e-10,
)
quats_bolted = sol_to_yup_quats(sol_bolted)
renderer_bolted, action_bolted = make_top_animation(quats_bolted, sol_bolted.t, show_bolt=True)
renderer_bolted
Principal Axes#
The inertia matrix \(\mathbf{M}_{b/c}\) is real and symmetric, so it has three orthogonal eigenvectors, the principal axes, and three positive eigenvalues, the principal moments of inertia. Spinning about a principal axis is special. The angular momentum \(\mathbf{h} = \mathbf{M}\boldsymbol{\omega}\) is then parallel to \(\boldsymbol{\omega}\), producing a steady spin with no wobble.
For the balanced top the spin axis \(z\) is already a principal axis (the inertia matrix is diagonal in the body frame). Adding the bolt breaks this. \(\mathbf{M}_\text{bolted}\) gains off-diagonal terms and its principal axes are tilted relative to the symmetry axis.
eigenvalues, eigenvectors = np.linalg.eigh(M_boltedtop)
print("Principal moments of inertia (kg·m²):")
print(np.round(eigenvalues * 1e6, 4), " (×10⁻⁶)")
print("\nPrincipal axes (columns of eigenvector matrix):")
print(np.round(eigenvectors, 4))
Principal moments of inertia (kg·m²):
[225.0022 265.1311 294.9333] (×10⁻⁶)
Principal axes (columns of eigenvector matrix):
[[-0.1552 -0.9879 0. ]
[ 0. 0. 1. ]
[-0.9879 0.1552 0. ]]