Gyroscopic forces in ODE

Conservation of angular momentum is a bit questionable in ODE.  If your rigid bodies have inertial tensors whose principal axes are all equal, everything is fine.  Otherwise, you will get into trouble.  A freely rotating body will tend to gain energy rather rapidly, eventually spinning fast enough to crash the simulator.  This was discussed somewhat recently on the ODE mailing list. I have generally worked around this problem by either adding global angular "friction" to the system or applying a maximum threshold to angular velocity.  This is fine but a little ugly.  You can also disable the gyroscopic term for each body which will result in conservation of angular velocity instead of angular momentum: a spinning figure skater pulling her arms in won't experience any acceleration.

The newsgroup discussion referenced a paper by Claude Lacoursière that addresses this problem (building on a paper by Anitescu and Potra).  The essence of the paper seems to be that you can achieve stability by manipulating the inertia/mass matrix to account for the gyroscopic torques.

It seems worthwhile to start with the computations affecting angular velocity in ODE's time step.  I'll use the notation in Lacoursière's paper.  Before resolving constraints for each time step of $h$ seconds, ODE first computes the inverse inertia matrix for time $n$ in the global frame: $\mathcal{I}^{-1}_{n} = R_n\mathcal{I}^{-1}_{0}R^{\intercal}_n$.  If gyroscope mode is enabled, it then computes the extra gyroscopic torques for the time step: $\tilde{\tau}_n = [{\mathcal{I}_n}\omega_n]\times\omega_n + \tau_n$.  Gravitational forces are also computed at this time, but they're not relevant for this discussion.

After computing the constraint Jacobian, ODE prepares a system for the LCP solver.  The solution is effectively found in "acceleration space".  Essentially, we start with $F=ma$ and then find the $F$ that satisfies $m^{-1}F=a$ where $a$ is the acceleration necessary to satisfy our constraints.  For rotation, the conserved value will be computed as $\mathcal{I}^{-1}_{n}\tilde{\tau}_n+\frac{\omega_n}{h}.$  Once the solver finds the constraint forces, the new angular velocity will be computed:

$\omega_{n+1} =\omega_{n}+h\mathcal{I}^{-1}_{n}\tilde{\tau}^*_n$ (where $\tilde{\tau}^*_n$ includes the constraint forces).
Ignoring the constraint forces, this expands to the explicit time step:
$\omega_{n+1} = \omega_n+h\mathcal{I}^{-1}_{n}[\hat{L}_n\omega_n+\tau_n]$ where $\hat{L}_n$ is the angular momentum, $\mathcal{I}_n\omega_n$, in "cross-product-matrix" form.  As mentioned before, this explicit form gains energy much like an explicit spring.  That is bad.

The paper goes on to compute the implicit form in body-local coordinates and finds that it's equivalent to solving the explicit form with a different inertia matrix: $\tilde{\mathcal{I'}}_n=\mathcal{I}_0-h\hat{L'}_n$.  According to the paper, this "is rather bizarre since then, the modified inertia tensor is non-symmetric" (pg. 8).  The paper shrugs at this form, stating that it would cause problems for most multibody frameworks and is therefore not very interesting, but they still give a function for the stepper:
$\omega_{n+1}'=[I+h\alpha C_n +h^2\alpha C^2_n]\omega_n' + h[I+h\alpha C_n + h^2\alpha C^2_n]\mathcal{I}^{-1}_0\tau_n'$
where
$\alpha = 1/(1+h^2\kappa)$
$\kappa = \det(\mathcal{I}^{-1}_0)\hat{L'}^{\intercal}_n\mathcal{I}_0\hat{L'}_n$
and
$C_n = \mathcal{I}^{-1}_0\hat{L'}_n$

I derived things a little differently.  ODE solves things in global space, rather than local space.  If we start with eq. 5.2 from the paper: $\tilde{\mathcal{I'}}_n\omega'_{n+1}=\mathcal{I}_0\omega'_n+h\tau'$ and move things into global space, we get the following:
$\omega_{n+1}=\tilde{\mathcal{I}}^{-1}_n\mathcal{I}_n\omega_n+h\tilde{\mathcal{I}}^{-1}_n\tau_n$
If we set the right-hand side to be equal to the explicit computation and then solve for $\tilde{\tau}_n$, we get $\tilde{\tau}_n=\frac{1}{h}\mathcal{I}_n[\tilde{\mathcal{I}}^{-1}_n\mathcal{I}_n\omega_n+h\tilde{\mathcal{I}}^{-1}_n\tau_n-\omega_n]$.
This is a lot uglier than $\tilde{\tau}_n = [{\mathcal{I}_n}\omega_n]\times\omega_n + \tau_n$, but as long as $\tilde{\mathcal{I}}_n$ is invertible, ODE shouldn't have much trouble with it.  The LCP solver never needs to see the funky mass matrix.  If we group terms a little better, we can see that there's not all that much extra to compute:
$\tilde{\tau}_n=\mathcal{I}_n\tilde{\mathcal{I}}^{-1}_n (\frac{1}{h}L_n)+\mathcal{I}_n\tilde{\mathcal{I}}^{-1}_n\tau_n-(\frac{1}{h}L_n)$
We already need to compute $\mathcal{I}_n=R_n\mathcal{I}_{0}R^{\intercal}_n$.
$L_n = \mathcal{I}_n\omega_n$ is cheap.
$\tilde{\mathcal{I}}_n=\mathcal{I}_n-h\hat{L}_n$ is cheaper.
I believe that inverting a 3x3 matrix can be done without too much cost and my feeling is that this matrix should always be invertible, although the matter merits additional attention.

I tested a straightforward implementation in ODE and found that it takes 3 to 7 times longer than the simple matrix-multiply and cross-product of the explicit method, but it's stable instead of exploding and that computation is still pretty small compared to the constraint solver or collision detection.

Local angular velocity of a single spinning body in ODE.  Inertia tensor is diag(10,4,2).  Initial angular velocity is (1,2,3).  Time step is h=1/60.  With gyroscope disabled, the body experiences conservation of angular velocity.  With explicit gyroscopic torques, it explodes.  With implicit computations, it converges.  It ought to maintain a closed fixed orbit that looks like the explicit one.

If we plot the local angular velocity of a spinning body, we get something that looks a lot like Fig. 8.1. from Lacoursière's paper (but cooler with animation).

Interestingly, the paper's justification for their inertia tensor of choice is

...it has three different values but is not an extreme case of a long thin object.  This illustrates many of the difficulties that can be encountered in practice. (pg. 11)

 However, when I went to visualize the geometry of an object with such an inertia tensor, I found that I couldn't do it.  It's impossible to create a box or ellipsoid of uniform density with that tensor unless you allow it to have imaginary dimensions.  My intuition says that you can't create any object with that tensor.  Stack-exchange verifies that the principal axes of an inertia tensor need to satisfy the triangle equality. 

    If, instead, we create a box with the density of water using the same ratio for dimensions: $[1,.4,.2]\text{m}$, then we get an $80\text{kg}$ box with inertia tensor $I_0\approx\text{diag}(1.33,6.93,7.73)$.  The end effect with this mass moment is pretty similar: the implicit method converges to a stable point and the explicit explodes.

I tried averaging the two methods or alternating between methods on successive passes, but, at least for this set of conditions, the trajectory was almost the same as the implicit trajectory.  The paper discusses other approaches that are "symplectic and time-reversible" and are much better at approximating the correct behavior, but these seems to use completely different integration methods and might not be amenable to implementation with popular rigid-body physics engines.

One remaining issue that should be addressed is the fact that I used a matrix inverse in my computation.  As I have experimented with it, I've found that there's a danger of encountering numeric instabilities with regular masses.