How does the orientation quaternion evolve in time under a measured angular velocity? This builds on the Hamilton-Convention (scalar-first $\mathbf q=[q_w,\mathbf q_v]$, orientation $\mathbf q=\mathbf q_{GL}=\mathbf q_{B\to I}$), and leads directly into the prediction step of the 3D-ESKF.
The orientation of a rigid body is a unit quaternion
$$\mathbf q = q_w + q_x\,\mathbf i + q_y\,\mathbf j + q_z\,\mathbf k = \begin{bmatrix}q_w\\\mathbf q_v\end{bmatrix}, \qquad \lVert\mathbf q\rVert = 1 .$$
The unit-norm constraint is essential — only unit quaternions represent rotations — and we will have to maintain it when integrating (below).
For a very small rotation angle $\Delta\phi$ (over a small interval $\Delta t$) about a unit axis $\mathbf u$, the incremental rotation quaternion is
$$\Delta\mathbf q_{\text{rot}} = \cos\!\Big(\tfrac{\Delta\phi}{2}\Big) + \mathbf u\,\sin\!\Big(\tfrac{\Delta\phi}{2}\Big).$$
For small angles, $\cos(\Delta\phi/2)\approx1$ and $\sin(\Delta\phi/2)\approx\Delta\phi/2$, so
$$\Delta\mathbf q_{\text{rot}} \approx 1 + \tfrac12\,\mathbf u\,\Delta\phi = \begin{bmatrix} 1 \\ \tfrac12\,\Delta\boldsymbol\phi \end{bmatrix}, \qquad \Delta\boldsymbol\phi = \mathbf u\,\Delta\phi\ \ (\text{rotation vector}).$$
Because the state is the local→global orientation $\mathbf q=\mathbf q_{GL}$, the axis $\mathbf u$ (and the rotation vector $\Delta\boldsymbol\phi$) are expressed in the local (body) frame.
The angular velocity is the rate of turning about the instantaneous axis, $$\vec\omega(t) = \frac{d\phi}{dt}\,\mathbf u(t) = \begin{bmatrix}\omega_x\\\omega_y\\\omega_z\end{bmatrix},$$ so over a small step $\vec\omega\,\Delta t = \mathbf u\,\Delta\phi = \Delta\boldsymbol\phi$. Written as a pure quaternion, $$\boldsymbol\omega = \begin{bmatrix}0\\\vec\omega\end{bmatrix},$$ the small-angle increment becomes $$\Delta\mathbf q_{\text{rot}} \approx 1 + \tfrac12\,\boldsymbol\omega\,\Delta t .$$ A gyroscope is fixed to the body, so the measured rate is in the body frame, $\vec\omega=\vec\omega_B$.
The orientation at $t$ is $\mathbf q(t)=\mathbf q_{GL}$. We update it by composing with the body-frame increment. In the Hamilton convention a rotation expressed in the local frame multiplies on the right:
$$\mathbf q(t+\Delta t) = \mathbf q(t)\otimes\Delta\mathbf q_{\text{rot}} .$$
Dual case (for reference). If instead the angular velocity were given in the world frame $\vec\omega_W$, the increment would multiply on the left: $\mathbf q(t+\Delta t)=\Delta\mathbf q_{\text{rot}}\otimes\mathbf q(t)$, giving $\dot{\mathbf q}=\tfrac12\boldsymbol\omega_W\otimes\mathbf q$. We use the body-frame (right) form throughout (cf. Solà §7 for the global-error variant).
Taking $\Delta t\to0$ in $\mathbf q(t+\Delta t)=\mathbf q(t)\otimes(1+\tfrac12\boldsymbol\omega\Delta t)$ gives the quaternion kinematics:
$$\boxed{\;\dot{\mathbf q} = \tfrac12\,\mathbf q\otimes\boldsymbol\omega\;}$$
Since the quaternion product is bilinear, this is also a linear ODE in $\mathbf q$, $$\dot{\mathbf q} = \tfrac12\,[\boldsymbol\omega]_R\,\mathbf q,\qquad [\boldsymbol\omega]_R = \begin{bmatrix} 0 & -\omega_x & -\omega_y & -\omega_z\\ \omega_x & 0 & \omega_z & -\omega_y\\ \omega_y & -\omega_z & 0 & \omega_x\\ \omega_z & \omega_y & -\omega_x & 0\end{bmatrix},$$ where $[\boldsymbol\omega]_R$ is the right product-matrix of the pure quaternion $\boldsymbol\omega$. This $4\times4$ skew form is what the ESKF linearises.
(a) First-order (Euler). Approximate $\dot{\mathbf q}$ over one step: $$\Delta\mathbf q_{\text o}\approx\dot{\mathbf q}\,\Delta t = \tfrac12\,\mathbf q_k\otimes\boldsymbol\omega\,\Delta t, \qquad \mathbf q_{k+1} = \mathbf q_k + \Delta\mathbf q_{\text o} = \mathbf q_k + \tfrac12\,\mathbf q_k\otimes\boldsymbol\omega\,\Delta t .$$ This is the additive update; it is only an approximation and pushes $\mathbf q$ off the unit sphere, so it must be renormalised each step ($\mathbf q_{k+1}\leftarrow\mathbf q_{k+1}/\lVert\mathbf q_{k+1}\rVert$).
(b) Exact (zeroth-order hold). If $\boldsymbol\omega$ is assumed constant over the step, the increment is the exact exponential, which is already unit-norm: $$\mathbf q_{k+1} = \mathbf q_k\otimes\mathrm{Exp}(\boldsymbol\omega\,\Delta t),\qquad \mathrm{Exp}(\boldsymbol\omega\Delta t) = \begin{bmatrix}\cos(\lVert\vec\omega\rVert\Delta t/2)\\[2pt] \dfrac{\vec\omega}{\lVert\vec\omega\rVert}\sin(\lVert\vec\omega\rVert\Delta t/2)\end{bmatrix}.$$ Form (b) needs no renormalisation and is more accurate for the same $\Delta t$; form (a) is the linearisation that underlies the EKF/ESKF prediction.
We integrate a constant body-rate $\boldsymbol\omega$ and check: (1) the ODE $\dot{\mathbf q}=\tfrac12\mathbf q\otimes\boldsymbol\omega$ against a finite difference, (2) the matrix form $\mathbf q\otimes\boldsymbol\omega=[\boldsymbol\omega]_R\mathbf q$, (3) that un-normalised Euler drifts off unit norm while the exact scheme stays on it, and (4) the small orientation error of Euler+renormalise vs. the exact integrator.
import numpy as np
def qmul(a, b):
aw, av = a[0], a[1:]; bw, bv = b[0], b[1:]
return np.concatenate([[aw*bw - av @ bv], aw*bv + bw*av + np.cross(av, bv)])
def qconj(q): return np.concatenate([[q[0]], -q[1:]])
def qnorm(q): return np.linalg.norm(q)
def qexp(phi): # Exp(rotation vector) -> unit quaternion
ang = np.linalg.norm(phi)
if ang < 1e-12: return np.array([1.0, 0.0, 0.0, 0.0])
return np.concatenate([[np.cos(ang/2)], (phi / ang) * np.sin(ang/2)])
def Omega_R(w): # right product-matrix of the pure quaternion [0,w]
wx, wy, wz = w
return np.array([[0,-wx,-wy,-wz],[wx,0,wz,-wy],[wy,-wz,0,wx],[wz,wy,-wx,0]])
q0 = qexp(np.array([0.2, -0.5, 0.7])) # initial orientation
w = np.array([0.3, 1.2, -0.6]) # constant body-frame angular velocity [rad/s]
wq = np.concatenate([[0.0], w]) # ω as a pure quaternion
# (1) ODE vs finite difference of the exact solution q(t) = q0 ⊗ Exp(ω t)
qt = lambda t: qmul(q0, qexp(w * t)); h = 1e-6
print("(1) ||q̇ - ½ q⊗ω|| :", np.linalg.norm((qt(h) - qt(-h)) / (2*h) - 0.5*qmul(q0, wq)))
# (2) matrix form
print("(2) ||q⊗ω - [ω]_R q|| :", np.linalg.norm(qmul(q0, wq) - Omega_R(w) @ q0))
# (3)-(4) integrate for 10 s
dt, T = 0.02, 10.0; n = int(T / dt)
q_euler = q0.copy(); q_norm = q0.copy(); q_exact = q0.copy()
for _ in range(n):
q_euler = q_euler + 0.5*qmul(q_euler, wq)*dt # Euler, no renormalisation
q_norm = q_norm + 0.5*qmul(q_norm, wq)*dt
q_norm /= qnorm(q_norm) # Euler + renormalise
q_exact = qmul(q_exact, qexp(w*dt)) # exact (constant ω)
print("(3) Euler(no renorm) ||q|| @10s:", round(qnorm(q_euler), 4), " | exact ||q|| =", round(qnorm(q_exact), 4))
dq = qmul(qconj(q_exact), q_norm)
print("(4) Euler+renorm vs exact @10s :", round(np.degrees(2*np.arctan2(np.linalg.norm(dq[1:]), abs(dq[0]))), 3), "deg")
Normalisation note. Because the first-order update $\mathbf q_{k+1}=\mathbf q_k+\tfrac12\mathbf q_k\otimes\boldsymbol\omega\Delta t$ is only an approximation, $\mathbf q$ slowly leaves unit length (the demo shows $\lVert\mathbf q\rVert\!\approx\!1.05$ after 10 s) — which corrupts the rotation it represents. Renormalise every step, or use the exact exponential update, which stays on the unit sphere by construction.
Even with exact integration and renormalisation, integrating gyro data alone ("dead reckoning") drifts, because we integrate the measured rate, which carries:
The remedy is sensor fusion: an accelerometer (gravity) and/or magnetometer provide an absolute orientation reference that corrects the drift, and a gyro bias state is estimated on-line. Combining the kinematics above (prediction) with these corrections is exactly the error-state Kalman filter.
Literature