Quaternion calculus and angular velocity

I’m having some trouble with my quaternion rotation class that I thought I’d figured out long ago, but now that I’m carefully checking the numbers, it’s not working as I expect. The issue is I want to move between rotations and angular velocities easily, and my rotation class uses quaternions for the internal representation. The math is straightforward and laid out in many different references:

w = ( 2 / dt * ( q1 - q0 ) * q0^-1 ).xyz;

q1 = q0 + dt * 1/2 * ( <wx,wy,wz,0> * q0 )

These equations relate an initial rotation quaternion q0, a final rotation quaterntion q1, a timestep dt, and an angular velocity 3-vector w. So now let’s look at an example in detail.

// identity rotation
q0 = <0,0,0,1>

// 90 degrees about x-axis
q1 = <0.707,0,0,0.707>

dt = 1

From the first equation (derivative, outputs angular velocity), I’d expect to get <1.571,0,0> - 90 degrees about the x-axis per unit time. But let’s step through the equation:

// component-wise subtraction
q1 - q0 = <0.707,0,0,-0.293>

// inverse of identity is identity
q0^-1 = q0

// multiply by identity
( q1 - q0 ) * q0^-1 = <0.707,0,0,-0.293>

2 / dt = 2

// multiply by scalar
2 * ( ( q1 - q0 ) * q0^-1 ) = <1.414,0,0,-0.586>

So w = <1.414,0,0>, which is 81 degrees about the x-axis! The same problem happens when integrating. If I set w = <1.571,0,0> (90 degrees about x-axis), dt = 1, and q0 = <0,0,0,1>, the output q1 = <0.618,0,0,0.786>, which is 76 degrees about the x-axis.

And of course, the sanity check is that I should be able to take the derivative and then the integral and get the original result back, but if I start with

q0 = <0,0,0,1>
q1 = <0.707,0,0,0.707>
dt = 1

I get after the derivative and then integral

w = <1.414,0,0>
q1’ = <0.577,0,0,0.816>

Which is only 70 degrees about the x-axis (so the errors from derivation and integration aren’t symmetric - they accumulate instead of cancelling out).

What’s going on? Anyone have some insight? Thanks!

-stephen diverdi

Hm, I believe I’ve figured out what the problem is. The math is fine, the issue is that the quaternion derivatives are only valid for small rotations (small timesteps and small angular velocities). If I change my example to use an angular velocity of 1 degree per unit time and integrate / derive over a timestep of 1, the math works out fine.

Conceptually this makes sense too. Only quaternions that are of unit length (on the surface of the hypersphere) represent rotations. Imagine this in 2D to make it easier - a valid rotation is a point on the unit circle. The angular velocity says how fast the point moves around the perimeter of the circle. However, the derivative gives the tangential velocity of the point - moving any amount in either direction along the tangent takes the point off the unit circle. For small movements and with renormalization, the error this causes is very small. As the size of the movement increases, the error increases significantly. Like so many other techniques, this approach just treats rotations as locally linear and then takes small steps.

I guess there’s probably a more accurate way to do this (integrate large rotations) without resorting to taking many small steps by moving the point (quaternion) along the hypersphere, rather than tangential to it. I wonder if the math for this is floating around somewhere?..

Thanks for looking,

-stephen diverdi