How to simulate rotational instability?

Question

I'm trying to simulate (for an educational game) the well-known effect that rotating objects with three nonequal moments of inertia are unstable when rotated around the middle axis.

Some explanations of this effect are given here:

• Stability of rotation of a rectangular prism
• Why does this object periodically turn itself?
• https://mathoverflow.net/questions/81960
• http://farside.ph.utexas.edu/teaching/336k/Newtonhtml/node71.html

Most of this math is, alas, over my head (though I'm trying). But for my purposes, I need more than just "is the form of a harmonic oscillator" or "is not a harmonic oscillator"; I need to actually reproduce the effect in simulation.

I attempted this using Unity's built-in physics engine here:

You can try it yourself if you have (or install) the free Unity plugin; it shows a 1x4x7 block of uniform density, rotating about its middle axis. With the "Poke" button you can induce a small random torque. Poking the block repeatedly can knock its axis askew -- but once you stop poking it, the axis stays put, and it rotates steadily around whatever direction that is. Under no circumstances have I been able to make it tumble (as seen in this video of a deck of cards, or this one of some simulation).

And its lack of tumbling makes perfect sense to me. As I understand it, the state of a rigid body is defined by its position, velocity, orientation, and angular velocity. Absent any external forces, the velocity and angular velocity should remain unchanged. Angular velocity can be described as an axis and a speed of rotation. So, without any external forces acting on it, how can the axis of rotation possibly change?

Clearly there's something missing both from my intuitive understanding, and from the physics engine within Unity. Don't focus too much on the latter; I can program my own physics engine, if I understand what it should do. What's the key bit I'm missing, that explains how (and in what manner) the axis of rotation can change without any external force? In pseudocode, simple forward Euler integration style, how would one simulate this?

Answer 1

$\vec\omega = I^{-1} \vec L$, and $\vec L$ is constant in the absence of external forces. The bit that I think you're missing is that $I$ rotates with the rigid body, so it is not constant in general and neither is $\vec\omega$.

I played with your online example, and the angular velocity does seem to always remain constant when I'm not poking the block, which is consistent with an inertia tensor that's a scalar multiple of the identity. I've never used Unity, but it seems to support arbitrary tensors of inertia via Rigidbody.inertiaTensor and Rigidbody.inertiaTensorRotation, so maybe you just need to set those properly.

If you end up having to roll your own physics, here's some naive, untested, and potentially wrong pseudocode:

const double time_step = ...;
const Quaternion L = {0, Lx, Ly, Lz};  // angular momentum
const double K1 = 1/I1, K2 = 1/I2, K3 = 1/I3;  // reciprocals of principal-axis moments of inertia; no idea if there's a standard letter for this
Quaternion orientation = {1, 0, 0, 0};
while (1) {
    Quaternion transformed_L = conjugate(orientation) * L * orientation;
    Quaternion transformed_omega = {0, K1 * transformed_L.i, K2 * transformed_L.j, K3 * transformed_L.k};
    Quaternion omega = orientation * transformed_omega * conjugate(orientation);
    orientation = quaternion_exp(time_step * omega) * orientation;
    // ... or orientation = normalize((1 + time_step * omega) * orientation);
    // ... or orientation = normalize((1 + time_step * omega + 0.5 * (time_step * omega)**2) * orientation);
    output(orientation);
}