# How to simulate rotational instability?

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:

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?

-
The angular momentum of an object doesn't change when an object tumbles, but the actual orientation of the object with respect to the axis the object is rotating around does. I believe you have identified a problem with Unity's built-in physics engine. – Peter Shor Jul 28 '14 at 17:59
You need your moment of inertia tensor to move with the object. For a rotating object, the instantaneous values of the moment of inertia along X, Y and Z change, and so the rotation about X, Y and Z must change in order to maintain constant angular momentum. This is why there is tumbling. – Floris Jul 28 '14 at 19:38
OK, this is almost making sense. If we think about the rotation as decomposed into individual rotations around the world X, Y, and Z axes, then yeah, the distance of the various particles in the body to these axes is changing, so (as you say) the rates of rotation would need to change to maintain angular momentum. For that matter, the same argument applies if we use the object's local coordinate system: once the rotation axis is perturbed, we'd have to recompute the moments of inertia relative to this new axis. Hmm. I don't yet understand it well enough to code it, but I'm getting closer! – Joe Strout Jul 29 '14 at 3:03
I've tried coding up my understanding here (forum.unity3d.com/threads/259514); I've also updated the web demo (luminaryapps.com/temp/RotationDemo). It seems to work, but I'm a little uncertain about my conversion of I_inverse from local to world coordinates. – Joe Strout Jul 29 '14 at 15:49

$\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);
}

-
Those are set automatically from the geometry (or supposed to be). I had already checked inertiaTensor, and it's (54.2, 41.7, 14.2), which matches expectations. I hadn't checked inertiaTensorRotation, but it's the identity quaternion (i.e. no rotation). – Joe Strout Jul 29 '14 at 2:51
I added additional logging of these at other times, and they never change. I'm not sure now whether that is expected or not. I'm still trying to digest your first comment above. – Joe Strout Jul 29 '14 at 2:53
inertiaTensor looks fine, and should be constant. The orientation of the tensor should be the same as the orientation of the object. I can't tell whether inertiaTensorRotation is relative to the canonical orientation or relative to the object's current orientation (as given by Rigidbody.rotation), In the former case it should be equal to Rigidbody.rotation. In the latter case it should be the identity quaternion. So I hope the former is correct since otherwise I have no idea what's wrong. – benrg Jul 29 '14 at 3:06
After more digging, I believe inertiaTensorRotation is relative to the object's current rotation, and that for it to be identity is normal. Though just for a lark, I tried setting it equal to the Rigidbody.rotation on each physics frame. This causes the physics engine to freak out a bit in other ways (actually translating the object periodically), but still doesn't cause any tumble in the rotation axis. But I think this is abusing the system -- I don't think it's intended for inertiaTensorRotation to change during the simulation. – Joe Strout Jul 29 '14 at 4:12
By the way, I've learned that Unity uses the Nvidia Physx engine, documented (thinly) here (docs.nvidia.com/gameworks/content/gameworkslibrary/physx/…). Instead of the usual 3x3 inertia tensor matrix, they keep just the diagonal of this matrix, and say: "If you have a non diagonal inertia tensor matrix, you need to diagonalize it and set an appropriate mass space transform." I believe this is what inertiaTensorRotation is for. Oh, and these are all relative to the object's coordinate system, so it's not surprising that they're constant. – Joe Strout Jul 29 '14 at 4:15