Mechanics 2: Linear and Rotational Dynamics - 3D Math Primer for Graphics and Game Development - page 638

(component-wise) derivative of the orientation, take a simple Euler step
on each component independently, and then correct the orientation. With
quaternions, the derivative is given by

Derivative of orientation
quaternion q for an
object with angular
velocity ω

q = 1

2 ω q ,
where the 3D angular velocity vector ω has been extended into quaternion
space with w = 0. (Eberly [17, Section 10.6] derives this result; we will
not prove it here.) Note that we do not expect q to be a rotation (unit)
quaternion. Neither do we expect the result of Euler integration, q k+1 ≈
q k + h q k to have unit length, and thus it must be normalized.
The technique of integration of orientation just described is a standard
one. When using this technique, there are two sources of error. The first
is caused by the Euler integration itself, in which we ignore angular ac-
celeration (and higher derivatives) and proceed as if the angular velocity
were constant. This error existed with linear data as well as angular; but
when integrating linear quantities it was the only source of error. The sec-
ond source of error is due to the use of component-wise derivatives, which
does not take into consideration the interdependence of the components of
the rotation matrix or quaternion. This type of error is unique to angular
data, since the components of positional data are independent (at least
when Cartesian coordinates are used). Fortunately, this source of error can
be eliminated.
Assume for the moment that the object whose orientation is described
by the rotation matrix R or quaternion q is rotating with constant angular
velocity ω. In this common case, no accuracy is lost by ignoring angular
acceleration, but there is a loss in accuracy, which increases with larger hω.
The solution is straightforward: determine the finite rotation that would oc-
cur in this time step, and then apply the appropriate angular displacement.
We already have the tools at our disposal. We convert the rotation hω from
exponential map form into axis-angle form (see Section 8.4). This angular
displacement can then be converted to a rotation matrix (see Section 5.1.3)
or quaternion (see Section 8.5.2), and concatenated with the current ro-
tation. Essentially, what we have done is to choose a better coordinate
system in which to perform the Euler integration.
Since it is slightly more expensive to perform the second method, a valid
question is: Is it worth it? In the common situation of constant angular
velocity, the second method is exact for any step size h, ignoring errors due
to floating-point roundoff. In this case, if the angular velocity is high, or
the accuracy required is high (e.g., the hands of a clock), then switching to
this alternative method will probably be a win. However, the error intro-
duced by Euler integration can interfere with the error introduced by the
component-wise derivative method either constructively or destructively,