Introduction

A while ago I wanted to know the moment of inertia of a tetrahedron. I'd forgotten some of the basic stuff, and the calculation was a bit fiddly so I thought I'd write it up on here. There are three related articles:

  1. Some basic results.1
  2. A toy problem: the sphere.2
  3. The final calculation3 (this article).

A first attack

Following the line we sketched when thinking about the sphere, we'll proceeed by integrating \(\textbf{r}\textbf{r}^T\) over the body:

\[ \textbf{C} = \int \;\ \textbf{r}\; \textbf{r}^T \; \mu \; dV. \]

This looks simple enough, but it's messy. In particular, if we treat it as three nested integrals over Cartesian axes, the limits of the integration will be fiddly to get right.

So, let's consider a simpler case.

A special case

Consider a tetrahedron with vertices \((0,0,0)\), \((1,0,0)\), \((0,1,0)\), and \((0,0,1)\) in Cartesian coordinates \(\textbf{u} = (u,v,w)\). Suppose further that it has unit uniform density.

A moment of thought will show that the volume inside the tetrahedron is defined by these four inequalities:

\[ \begin{align} u & \ge 0, \\ v & \ge 0, \\ w & \ge 0, \\ u + v + w & \le 1. \end{align} \]

Accordingly the integral over the tetrahedron is just

\[ \int_0^1 \; \int_0^{1-u} \; \int_0^{1 - u - v} dw \; dv \; du \]

It's worth noting that despite appearances this is symmetric in \(u\), \(v\), and \(w\).

In the following calculations we'll use subscript \(s\) to indicate that we're talking about the special case.

Volume

The integral above is useful: it's just the volume of the tetrahedron:

\[ \begin{align} V_s &= \int_0^1 \; \int_0^{1-u} \; \int_0^{1 - u - v} dw \; dv \; du, \\ &= \int_0^1 \; \int_0^{1-u} (1 - u - v) dv \; du, \\ &= \int_0^1 \; \frac{1}{2} (1 - u)^2 \; du, \\ &= \frac{1}{6}. \end{align} \]

Given that it has unit density, that's also its mass:

\[ M_s = \frac{1}{6}. \]

Centre of mass

Although we don't actually need this, let's calculate it anyway.

Given unit density, the centre of mass \(\textbf{p}_s\) is just:

\[ \textbf{p}_s = \frac{1}{M_s} \; \int \;\ \textbf{u}\; du \; dv \; dz. \]

We'll find that integral useful, so define

\[ \textbf{q}_s = \int \;\ \textbf{u}\; du \; dv \; dz, \]

Although it might seem that we have to do three integrals, one for each component of \(\textbf{u}\), recall that the tetrahedron is symmetric in \(u\), \(v\) and \(w\), so the answers will be the same. Let's just do the integral of \(u\), because it's the easiest:

\[ \begin{align} q_u = q_v = q_w &= \int_0^1 \; \int_0^{1-u} \; \int_0^{1 - u - v} u \; dw \; dv \; du, \\ &= \frac{1}{2} \int_0^1 \; u (1 - u)^2 \; du, \\ &= \frac{1}{24}, \\ \textbf{q}_s &= \frac{1}{24} \left(\begin{array}{c}1 \\ 1 \\ 1\end{array}\right). \end{align} \]

Accordingly,

\[ \begin{align} \textbf{p}_s &= \frac{1}{M_s} \textbf{q}_s, \\ &= \frac{1}{4} \left(\begin{array}{c}1 \\ 1 \\ 1\end{array}\right). \end{align} \]

Moment of Inertia

Let's start with

\[ \begin{align} \textbf{C}_s &= \int \textbf{u} \textbf{u}^T \; du \; dv \; dw \\ &= \int \left(\begin{array}{ccc} u^2 & uv & wu \\ uv & v^2 & vw \\ wu & vw & w^2 \end{array} \right) \; du \; dv \; dw \end{align} \]

and again note that symmetry dictates that there are only two distinct terms:

\[ \begin{align} C_{uu} = C_{vv} = C_{ww} &= \int_0^1 \; \int_0^{1-u} \; \int_0^{1 - u - v} u^2 \; dw \; dv \; du, \\ &= \int_0^1 \; \int_0^{1-u} u^2 \; (1 - u - v) \; dv \; du, \\ &= \frac{1}{2} \int_0^1 \; u^2 (1 - u)^2 \; du, \\ &= \frac{1}{60}. \end{align} \]

and

\[ \begin{align} C_{uv} = C_{vw} = C_{wu} &= \; \int_0^1 \; \int_0^{1-u} \; \int_0^{1 - u - v} u v \; dw \; dv \; du, \\ &= \int_0^1 \; \int_0^{1-u} (1 - u - v) \; u v \; dv \; du, \\ &= \frac{1}{6}\int_0^1 \; u (1 - u)^3 \; du, \\ &= \frac{1}{120}. \end{align} \]

So,

\[ \begin{align} \textbf{C}_s &= \frac{1}{120} \left(\begin{array}{ccc} 2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right), \\ &= \frac{M_s}{20} \left(\begin{array}{ccc} 2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right), \end{align} \]

and

\[ \begin{align} \textbf{I}_s &= \left(\textrm{Tr}\;\textbf{C}_s\right)\, \boldsymbol{\mathbb{I}} - \textbf{C}_s, \\ &= \frac{M_s}{20} \left(\begin{array}{ccc} 4 & -1 & -1 \\ -1 & 4 & -1 \\ -1 & -1 & 4 \end{array} \right). \end{align} \]

The general tetrahedron

Consider a general tetrahedron with vertices \(\textbf{a}, \textbf{b}, \textbf{c}, \textbf{d}\). We obviously could proceed as above, but the limits on the integrals will be messy.

Instead, consider this transformation:

\[ \begin{align} \textbf{x} &= \textbf{R} \textbf{u} + \textbf{a}, \\ &= \left(\begin{array}{c} \textbf{b} - \textbf{a} & \textbf{c} - \textbf{a} & \textbf{d} - \textbf{a} \end{array} \right) \textbf{u} + \textbf{a}. \end{align} \]

In other words the columns of \(\textbf{R}\) are just e.g. \(\textbf{b} - \textbf{a}\).

Although we've singled out \(\textbf{a}\) here, it's important to remember that the answer must be symmetric in all the vertices: after all how we label them is entirely arbitrary.

In particular consider these cases:

\[ \textbf{u} = \left(\begin{array}{ccc}0 \\ 0 \\ 0\end{array}\right), \, \left(\begin{array}{ccc}1 \\ 0 \\ 0\end{array}\right), \, \left(\begin{array}{ccc}0 \\ 1 \\ 0\end{array}\right), \, \left(\begin{array}{ccc}0 \\ 0 \\ 1\end{array}\right). \]

and convince yourself that these points transform to:

\[ \textbf{x} = \textbf{a}, \textbf{b}, \textbf{c}, \textbf{d}. \]

So, instead of integrating \(\textbf{x}\) over the general tetrahedron, we can simply integrate \(\textbf{u}\) over the particular case we considered above.

Of course, because there's a change of variable we'll need an appropriate Jacobian:

\[ dx\; dy\; dz = \left| \textbf{R} \right| du \; dv \; dz. \]

where \(\left| \textbf{R} \right|\) is the magnitude of the determinant of \(\textbf{R}\).

We can now calculate the volume, centre of mass, and moment of inertia for a general tetrahedron with just a bit of matrix algebra. In every case, the transformation above will transform the general integral over the tetrahedron to the special case we've already solved.

Volume

\[ \begin{align} V &= \int dx\; dy\; dz, \\ &= \left| \textbf{R} \right| \int du \; dv \; dz, \\ &= \frac{1}{6} \left| \textbf{R} \right|. \end{align} \]

A perfectly reasonable result when you consider that \(\left| \textbf{R} \right|\) is the triple scalar product of three sides of the tetrahedron.

We also note the the mass

\[ M = \frac{1}{6} \mu \left| \textbf{R} \right|. \]

Centre of mass

\[ \begin{align} \textbf{p} &= \frac{1}{M} \; \int \;\ \textbf{x}\; \mu \; dx \; dy \; dz, \\ &= \frac{\mu}{M} \left| \textbf{R} \right| \int \; \left( \textbf{R} \; \textbf{u} + \textbf{a} \right) \; du \; dv \; dw, \\ &= \frac{\mu}{M} \left| \textbf{R} \right| \left( \textbf{R} \int \; \textbf{u} \; du \; dv \; dw + \textbf{a} \int \; du \; dv \; dw \right) \\ &= 6 \left( \textbf{R} \textbf{q}_s + \textbf{a} V_s \right), \\ &= 6 \left( \frac{1}{24} \left(\textbf{b} + \textbf{c} + \textbf{d} - 3 \textbf{a}\right) + \frac{1}{6} \textbf{a} \right), \\ &= \frac{1}{4} \left(\textbf{a} + \textbf{b} + \textbf{c} + \textbf{d}\right), \\ &= \frac{1}{4} \sum_{i} \textbf{a}_i. \end{align} \]

Here the final sum is over the four vertices of the tetrahedron.

Moment of inertia

\[ \begin{align} \textbf{C} &= \int \; \textbf{x} \textbf{x}^T \; \mu \; dx \; dy \; dz, \\ &= \mu \left| \textbf{R} \right| \int \; \left( \textbf{R} \; \textbf{u} + \textbf{a} \right) \left( \textbf{R} \; \textbf{u} + \textbf{a} \right)^T \; du \; dv \; dw. \end{align} \]

If we expand the integrand there are effectively three terms to consider:

The term in \(\textbf{u} \textbf{u}^T\)

\[ \begin{align} \int \textbf{R} \textbf{u} \textbf{u}^T \textbf{R}^T \; du \; dv \; dw &= \textbf{R} \left( \int \textbf{u} \textbf{u}^T \; du \; dv \; dw \right) \textbf{R}^T \\ &= \textbf{R} \; \textbf{C}_s \; \textbf{R}^T, \\ &= \frac{1}{120}\left( 12 \textbf{a}\textbf{a}^T - 4\left(\textbf{a} \textbf{e}^T + \textbf{e} \textbf{a}^T\right) + \left(\textbf{b}\textbf{b}^T + \textbf{c}\textbf{c}^T + \textbf{d}\textbf{d}^T\right) + \textbf{e} \textbf{e}^T \right). \end{align} \]

where \(\textbf{e} = \left(\textbf{b} + \textbf{c} + \textbf{d}\right)\).

The term in \(\textbf{u}\)

\[ \begin{align} \int \textbf{R} \textbf{u} \textbf{a}^T \; du \; dv \; dw &= \textbf{R} \left( \int \textbf{u} \; du \; dv \; dw \right) \textbf{a}^T \\ &= \textbf{R} \; \textbf{q}_s \; \textbf{a}^T \\ &= \frac{1}{24} \left(\textbf{b} + \textbf{c} + \textbf{d} - 3 \textbf{a}\right) \textbf{a}^T. \end{align} \]

Of course there's the transpose of this term too.

The constant term

\[ \begin{align} \int \textbf{a} \textbf{a}^T \; du \; dv \; dw &= \textbf{a} \left( \int \; du \; dv \; dw \right) \textbf{a}^T \\ &= \textbf{a} \; V_s \; \textbf{a}^T \\ &= \frac{1}{6} \textbf{a} \textbf{a}^T. \end{align} \]

The whole integral

Putting these all together and collecting terms, gives this pleasing result:

\[ \begin{align} \int \; \left( \textbf{R} \; \textbf{u} + \textbf{a} \right) \left( \textbf{R} \; \textbf{u} + \textbf{a} \right)^T \; du \; dv \; dw &= \frac{1}{120} \left(\sum_i \textbf{a}_i \textbf{a}_i^T + \sum_{i} \textbf{a}_i \; \sum_{i} \textbf{a}_i^T \right), \\ \textbf{C} &= \frac{M}{20} \left(\sum_i \textbf{a}_i \textbf{a}_i^T + \sum_{i} \textbf{a}_i \; \sum_{i} \textbf{a}_i^T \right). \end{align} \]

So the trace is just

\[ \textrm{Tr}\;\textbf{C} = \frac{M}{20} \left(\sum_i \textbf{a}_i^T \textbf{a}_i + \sum_{i} \textbf{a}_i \; \sum_{i} \textbf{a}_i^T \right), \]

and thus

\[ \textbf{I} = \frac{M}{20} \left( \left(\sum_i \textbf{a}_i^T \textbf{a}_i + \sum_{i} \textbf{a}_i^T \; \sum_{i} \textbf{a}_i \right) \mathbb{I} - \left(\sum_i \textbf{a}_i \textbf{a}_i^T + \sum_{i} \textbf{a}_i \; \sum_{i} \textbf{a}_i^T \right) \right). \]

A more thoughtful solution

Although our approach above saves us from messy integrals, we still have to do some messy matrix algebra. We can do better!

It's clear from above that \(\textbf{C}\) is the sum of terms like \(\textbf{a}\textbf{b}^T\).

However we've noted before that the final expression for the moment of inertia is symmetric with respect to permutation of the vertices. That is, if we swap say \(\textbf{a}\) and \(\textbf{c}\), then the answer must not change.

So, the most general expression for \(\textbf{C}\) is:

\[ \textbf{C} = M \left( \alpha \sum_{i} \textbf{a}_i \textbf{a}_i^T + \beta \sum_{i} \textbf{a}_i \sum_{i} \textbf{a}_i^T \right). \]

where \(\alpha\) and \(\beta\) are unknown constants independent of the vertices.

Now, consider an infinitessimal tetrahedron, where all the vertices are the same, \(\textbf{h}\), say. For such a tetrahedron, it's clear that

\[ \textbf{C} = M \textbf{h} \textbf{h}^T \]

and so

\[ \begin{align} \textbf{h} \textbf{h}^T &= \left( \alpha \sum_{i} \textbf{h} \textbf{h}^T + \beta \sum_{i} \textbf{h} \sum_{i} \textbf{h}^T \right), \\ &= \left(4 \alpha + 16 \beta \right) \textbf{h} \textbf{h}^T \\ \beta &= \frac{1}{16} \left(1 - 4 \alpha \right). \end{align} \]

Substituting back into our expression for \(\textbf{C}\):

\[ \begin{align} \textbf{C} &= M \alpha' \left( \frac{1}{4} \sum_{i} \textbf{a}_i \textbf{a}_i^T - \left( \frac{1}{4} \sum_{i} \textbf{a}_i \right)\; \left( \frac{1}{4} \sum_{i} \textbf{a}_i^T \right) \right) + M \left( \frac{1}{4} \sum_{i} \textbf{a}_i \right)\; \left( \frac{1}{4} \sum_{i} \textbf{a}_i^T \right), \\ &= M \alpha' \left( \frac{1}{4} \sum_{i} \textbf{a}_i \textbf{a}_i^T - \textbf{p} \textbf{p}^T \right) + M \textbf{p} \textbf{p}^T. \end{align} \]

where \(\textbf{p}\) is the centre of mass: \(\frac{1}{4}\sum_i \textbf{a}_i\), and \(\alpha'\) is just a rescaled \(\alpha\).

It's clear that the two terms correspond to the covariance of the mass distribution, and the centre of mass itself—the latter has no \(\alpha'\) dependence.

However, this expression applies to any object with four equivalent vertices. It might be a solid tetrahedron, four isolated point masses, a hollow tetrahedral shell or so on. Each different type will have its own value of \(\alpha'\), but for a particular type we need only one scalar to calculate its moment of inertia.

For the solid tetrahedron, consider the special case from above, with vertices \((0,0,0)\), \((1,0,0)\), \((0,1,0)\), and \((0,0,1)\):

\[ \begin{align} \sum_{i} \textbf{a}_i \textbf{a}_i^T &= \left(\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right), \\ \textbf{p} \textbf{p}^T &= \frac{1}{16} \left(\begin{array}{ccc} 1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{array} \right). \end{align} \]

But we've already calculated \(\textbf{C}\) for this tetrahedron:

\[ \textbf{C}_s = \frac{M}{20} \left(\begin{array}{ccc} 2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right). \]

Thus for the solid tetrahedron

\[ \alpha' = \frac{1}{5}, \]

And so

\[ \begin{align} \textbf{C} &= \frac{M}{5} \left( \frac{1}{4} \sum_{i} \textbf{a}_i \textbf{a}_i^T - \textbf{p} \textbf{p}^T \right) + M \textbf{p} \textbf{p}^T, \\ &= \frac{M}{20} \left(\sum_i \textbf{a}_i \textbf{a}_i^T + \sum_{i} \textbf{a}_i \; \sum_{i} \textbf{a}_i^T \right). \end{align} \]

This cunning approach doesn't completely eliminate the integration, but we'd only have to integrate one component of \(\textbf{I}\) over the carefully chosen tetrahedron: the rest is really just symmetry.