Using optimization to find Principal axes of rigid body¶

My friend Rhett Allain posted an awesome video this week about numerically finding principal axes of a rigid body. His situation had masses confined to a plane and I wondered if I could use a full 3D approach for masses randomly scattered.

I'm going to do this two ways.

I'll keep rotating the points around and check to see if, in the new orientation, the moment of inertia tensor has zero off-diagonal elements. If so, then the normal axes after that rotation will be the principal axes. I can then watch to see where those axes go as I undo the rotation back to the original orientation.
I'll calculate the moment of inertia tensor in the original orientation (this'll have non-zero off-diagonal elements) and then look for a vector, $\vec{\omega}$ such that $\overleftrightarrow I\cdot\vec{\omega}$ yields a vector that's parallel to $\vec{\omega}$. That's really saying that the angular velocity vector and the angular momentum are parallel and that only happens along one of the principal axes. Then I'll look for another such situation in the plane that's perpendicular to the first one, and then I'll find the third axis by crossing the first two.

Both are an extension of what Rhett did, though I think (2) is a more direct analogy. Interestingly (2) is much less computationally intensive because (1) requires a constant recalculation of the intertia tensor (which means a bunch of sums) for every trial orientation of the system.

Method 1: look for an orientation that diagonalizes the inertia tensor¶

First I'll make some random points and shift them so the origin is at the center of mass:

In [1]:

pts=RandomReal[{-1,1},{10,3}];
com=Total[pts]/10;
pts=#-com&/@pts;
Show[Graphics3D[Point[pts]]]

Out[1]:

Next I'll calculate the three off-diagonal elements for a given set of points. The pattern after the colon basically says "hey this function should only evaluate if you give me a nested list of numbers" so that the minimizer below doesn't try to make it run symbolically.

In [5]:

Clear[offDiagonals]
offDiagonals[pts : { {__?NumberQ}.. }]:=Total[Table[pt[[#[[1]]]] pt[[#[[2]]]], {pt, pts}]]&/@{ {1,2}, {2,3}, {3,1}}

In [7]:

offDiagonals[pts]

Out[7]:

{-0.535606, 1.65996, 0.90678}

See, they're non-zero!

Next I want a formula for the rotation matrix for a given set of Euler angles ($\theta$, $\phi$, and $\psi$):

In [9]:

rot[theta_Real, phi_Real, psi_Real]:=RotationMatrix[phi, {0,0,1}].
                                     RotationMatrix[theta, {0,1,0}].
                                     RotationMatrix[psi, {0,0,1}];
rot[1.,2.,3.]//MatrixForm

Out[9]:

Ok, now lets see if NMinimize can find a set of Euler angles that does the trick.

In [12]:

{bestValue, EulerBest}=NMinimize[(tmp=offDiagonals[rot[theta,phi,psi].#&/@pts]).tmp,{theta,phi,psi}]

Out[12]:

Awesome! That took about 4 seconds to do and now we've got the angles we want. Here's an animation showing the rotation from the original orientation to the one that puts the principal axes along x, y, and z:

In [26]:

frame[frac_]:=(tmpPoints=rot[Sequence@@(frac {theta,phi,psi}/.EulerBest)].#&/@pts)
ListAnimate[Table[Show[Graphics3D[Point[frame[frac]],PlotRange->Sqrt[3] 1.1]], {frac,0.0,1.0,0.01}]]

Out[26]:

Now to get the axes in the original position, we just have to find where the cardinal axes go to if we invert the rotation (which, for these rotations, is just the transpose)

In [30]:

inverseRotation=Inverse[rot[Sequence@@({theta,phi,psi}/.EulerBest)]];
inverseRotationT=Transpose[rot[Sequence@@({theta,phi,psi}/.EulerBest)]];
inverseRotation//MatrixForm
inverseRotationT//MatrixForm

Out[30]:

See, I told you inverse and transpose did the same thing. Ok, now the correct axes:

In [42]:

newAxes=(inverseRotation.#&/@{ {1,0,0}, {0,1,0}, {0,0,1}});
newAxes//TableForm

Out[42]:

Here's a plot of the original points with those axes drawn in:

In [44]:

Show[Graphics3D[{Point[pts],
                 Red,
                 Arrow[{ {0,0,0}, #}]&/@newAxes}]]

Out[44]:

Method 2: find an $\vec{\omega}$ that's parallel with the angular momentum¶

We'll use those same initial points, but we'll start by calculating the full inertia tensor:

In [46]:

moi=Table[Total[Table[KroneckerDelta[i,j] pt.pt-pt[[i]] pt[[j]], {pt,pts}]], {i,3},{j,3}];
moi//MatrixForm

Out[46]:

Now we'll hunt for a vector that, when dotted into that matrix, produces a vector in the same direction.

In [49]:

{min, firstAxis}=NMinimize[(tmp=Cross[(moi.(omega={Sin[theta] Cos[phi], Sin[theta] Sin[phi], Cos[theta]})), omega]).tmp, {theta,phi,psi} ]

Out[49]:

Whoa! That was fast!

Ok, now we can rotate all the points to bring that vector to, say, the vertical and then look for another one that's in the x-y plane:

In [55]:

bringToVertical=RotationMatrix[{ {Sin[theta] Cos[phi], Sin[theta] Sin[phi], Cos[theta]}/.firstAxis, {0,0,1}}];
nextPts=bringToVertical.#&/@pts;
moi2=Table[Total[Table[KroneckerDelta[i,j] pt.pt-pt[[i]] pt[[j]], {pt,nextPts}]], {i,3},{j,3}];
{min, secondAxis}=NMinimize[(tmp=Cross[(moi2.(omega={Cos[phi], Sin[phi], 0})), omega]).tmp, {phi} ]

Out[55]:

Excellent! Now we can rotate that over to the x-axis and whatever's on the y-axis will be our third one. Then we can undo all the rotations and get the original axes:

In [60]:

bringToX=RotationMatrix[{ {Cos[phi], Sin[phi], 0}/.secondAxis, {1,0,0}}];
method2Axes=Inverse[bringToVertical].Inverse[bringToX].#&/@{ {1,0,0}, {0,1,0}, {0,0,1}}

Out[60]:

{{-0.352403, -0.797071, 0.490399}, {0.889413, -0.448265, -0.089453}, 
 
>   {0.291129, 0.404644, 0.866895}}

In [63]:

Show[Graphics3D[{Point[pts],
                 Red,
                 Arrow[{ {0,0,0}, #}]&/@method2Axes,
                 Blue,
                 Arrow[{ {0,0,0}, #}]&/@newAxes}]]

Out[63]:

Where I've put the method 2 ones in red and the method 1 ones in blue. They agree up to a random permutation.

Finally, let's compare with the Eigenvector approach:

In [67]:

Eigenvectors[moi]//TableForm

Out[67]:

In [68]:

newAxes//TableForm

Out[68]:

In [69]:

method2Axes//TableForm

Out[69]:

Yep, they all agree!

Your thoughts?¶

I'd love to hear what you think about this. Here are some starters for you:

This is cool! I really like . . .
This is dumb! If you could just do "Eigenvectors" all along why do all the rest?
Why do you never type two consecutive curly brackets?
How slow does method one get if you go up to 100 points? 1000?
Am I supposed to read all the axes horizontally or vertically?
So basically if Rhett posts something you're compelled to give it a try?
I'm in linear algebra right now and have no idea what eigenvectors are good for. Is this supposed to be helping me?