Use implicit Euler for gyroscopic torque

[Jondolf]

Objective

Fixes #400.

Gyroscopic motion is currently accounted for with a gyroscopic term computed using semi-implicit Euler. However, this sometimes leads to the torque just blowing up and making angular velocities essentially infinite, because semi-implicit Euler extrapolates velocities, and the gyroscopic torque is quadratic in angular velocity.

This manifested as angular velocities blowing up and rotation becoming invalid for objects with non-uniform inertia, especially when spinning at high speeds, as seen in the videos below:

2024-07-13.02-24-29.mp4

2024-07-13.02-23-23.mp4

Solution

Use implicit Euler for computing the gyroscopic torque in a much more stable way. The Newton Raphson method is used to solve the resulting non-linear equations. See Erin Catto's GDC 2015 slides on Numerical Methods for more information.

This fixes the big explosions and produces higher quality simulation:

2024-07-13.02-19-21.mp4

2024-07-13.02-30-28.mp4

Implicit Euler is slightly more expensive for this, but it should be acceptable in this case. We could explore other options in the future, like making gyroscopic motion optional, or trying an adaptive approach that only uses implicit Euler as a fallback.