Draft:SPIRAL algorithm

SPIRAL is a third-order integration algorithm for solving Euler's rigid body equations of motion using quaternions. It stands for Stable Particle Rotation Integration Algorithm. SPIRAL provides several advantages over existing methods, including improved stability, accuracy, and reduced overall computational cost. Furthermore, its formulation preserves the norm of the quaternion, eliminating the need for recurrent normalization.

This algorithm has practical applications in various fields, such as particle simulation, molecular dynamics (MD), discrete element methods (DEM), physics game engines, computer graphics, robotics, sensors, navigation systems, autonomous vehicles, and control systems.

Euler rigid body equations of motion
A rotating solid object is described by Euler's rigid body equations of motion:

$$\dot{\omega}_x = \frac{M_x}{I_x}+ \omega_y\,\omega_z \frac{I_y - I_z}{I_x}$$,

$$\dot{\omega}_y = \frac{M_y}{I_y}+ \omega_z\,\omega_x \frac{I_z - I_x}{I_y}$$,

$$\dot{\omega}_z = \frac{M_z}{I_z}+ \omega_x\,\omega_y \frac{I_x - I_y}{I_z}$$,

where $$\vec{\omega}(t)$$ is the angular velocity, and $$\vec{M}(t)$$ is the torque. Both quantities are in the body's principal axis reference frame. The body's principal moments of inertia are $$I_x$$, $$I_y$$, and $$I_z$$. Solving this system of non-linear first-order differential equations yields the body's angular velocity. Using the result, one can calculate the orientation of the body using the norm-conserving quaternion derivative:

$$\dot{q} = \frac{dq}{dt} = \frac{1}{2}\omega(t) q$$.

where $$\omega=\{0,\omega_x,\omega_y,\omega_z\} $$ is a pure imaginary quaternion representing the angular velocity $$\vec\omega=\{\omega_x,\omega_y,\omega_z\} $$ in the reference frame of the rotating body.

Algorithm
SPIRAL applies to both leapfrog and non-leapfrog architectures, so it's easily adaptable to different code bases. In this context, the leapfrog version.

First, using the known torques, update the angular velocity:

$$\vec{\omega}\left(t + \frac{\Delta t}{2}\right) = \vec{\omega}\left(t-\frac{\Delta t}{2}\right) + \frac{1}{6}\left(\vec{K}_1 + \vec{K}_2 + 4\,\vec{K}_3\right)$$,

with

$$\vec{K}_1 = \Delta t\,\dot{\vec{\omega}}\left(\vec{\omega}, \vec{M}(t)\right) $$.

$$\vec{K}_2 = \Delta t\,\dot{\vec{\omega}}\left(\vec{\omega} +\vec{K}_1, \vec{M}(t)\right) $$.

$$\vec{K}_3 = \Delta t\,\dot{\vec{\omega}}\left(\vec{\omega} + \frac{1}{4}\left(\vec{K}_1 + \vec{K}_2\right), \vec{M}(t) \right) $$.

The developers of SPIRAL proposed Strong Stability Preserving Runge-Kutta 3 (SSPRK3) to update angular velocity. Nevertheless, other algorithms can be utilized instead. It is important to note that for optimal performance, the integration algorithm selected must be at least as accurate as the quaternion one, making sure that the full advantage of SPIRAL's precision and stability can be taken. Also, as in other leapfrog algorithms, the angular velocity must be initialized half a step backwards. Note that torque is not recalculated for each Runge-Kutta stage, meaning that SPIRAL only requires one force calculation per time step.

Second, update the orientation:

$$q(t+\Delta t) = q(t)\cdot\left(\cos{\theta_1} + \frac{\vec{\omega}\left(t + \frac{\Delta t}{2}\right)}{\left|\vec{\omega}\left(t + \frac{\Delta t}{2}\right)\right|}\sin{\theta_1} \right),\quad \text{with} \quad \theta_1 \equiv \frac{\Delta t}{2}\left|\vec{\omega}\left(t + \frac{\Delta t}{2}\right)\right|\,. $$

Here the notation $$(a + \vec{b}) $$ represents a quaternion.

Bellow, a sample implementation in Julia of a time step:

Derivation
To derive an expression for $$q(t + \Delta t) $$, we start with the quaternion's time derivative:

$$\dot{q} = \frac{dq}{dt} = \frac{1}{2}\omega(t) q$$.

$$\omega=\{0,\omega_x,\omega_y,\omega_z\} $$ is a pure imaginary quaternion representing the angular velocity $$\vec\omega=\{\omega_x,\omega_y,\omega_z\} $$ in the reference frame of the rotating body. Assuming that it depends only explicitly on time and solve by separation of variables:

$$\int_t^{t + \Delta t} \frac{\text{d}q}{q} = \frac{1}{2}\int_{t}^{t + \Delta t} \omega(t^\prime) \text{d}t^\prime $$,

to obtain

$$q(t + \Delta t) = q(t)\,\exp\left(\frac{1}{2}\int_{t}^{t + \Delta t} \omega(t^\prime)\, \text{d}t^\prime\right) $$.

Expanding the angular velocity in a Taylor series around $$a $$ we get

$$\omega(t^\prime) = \omega(a) + \dot{\omega}(a)(t^\prime - a) + \frac{1}{2}\ddot{\omega}(a)(t^\prime - a)^2 + \mathcal{O}\left[\left(t^\prime - a\right)^3\right] $$.

Taking $$a=t + \frac{\Delta t}{2} $$ (leapfrog assumption) and perform the integral, we get

$$q(t + \Delta t) = q(t)\,\exp\left[\frac{\Delta t}{2}\omega\left(t + \frac{\Delta t}{2}\right) + \mathcal{O}(\Delta t^3)\right] $$,

with

$$O(\Delta t^3)\approx \frac{dt^3}{48}\ddot{\omega}\left(t + \frac{\Delta dt}{2}\right) $$.

The exponential of a purely imaginary quaternion $$\theta \hat{u} $$, where $$\theta $$ is a scalar and $$\hat{u}=\{0,u_x,u_y,u_z\} $$, a unit quaternion, reads

$$e^{\theta\hat{u}} = \cos{\theta} + \vec{u}\sin{\theta} $$,with the unit vector $$\vec u=\{u_x,u_y,u_z\} $$. Then

$$q(t+\Delta t) = q(t)\left(\cos{\theta_1} + \frac{\vec{\omega}\left(t + \frac{\Delta t}{2}\right)}{\left|\vec{\omega}\left(t + \frac{\Delta t}{2}\right)\right|}\sin{\theta_1} \right),\quad \text{with} \quad \theta_1 \equiv \frac{\Delta t}{2}\left|\vec{\omega}\left(t + \frac{\Delta t}{2}\right)\right| $$

Completing the derivation of the algorithm. This formulation preserves the norm of the initial quaternion, thus eliminating the need for subsequent normalization. Nonetheless, normalizing the quaternion every 50000 steps or so could prevent floating point precision-induced instabilities.

Accuracy and Stability
To evaluate the accuracy and stability of SPIRAL, its output can be compared against an analytical solution, providing a reliable means of assessing its performance relative to other integration algorithms. The chosen analytical solution serves as a benchmark, representing the rotational motion of a cylinder under the influence of a constant torque applied along one of its principal axes. The algorithms in the comparison are the algorithms used to integrate the rotational motion of particles in various well-established open-source and commercial discrete element method (DEM) programs like Yade, MercuryDPM , etc. The in the comparison are Runge-Kutta4 (but this one requires four force calculations per time step while the others only require one), Direct Euler, Velocity Verlet, Buss algorithm, Johnson algorithm , Fincham Leapfrog , Omelyan advanced leapfrog , algorithm found in MercuryDPM source code and both SPIRAL versions.First, the relative error between the analytical solution and the different integration algorithms as a function of the time step shows heuristically the convergence and accuracy of each algorithm. It is clear from this plot that all the algorithms, except for Runge-Kutta4 and SPIRAL, are second-order. SPIRAL is third-order, and Runge-Kutta4 is fourth-order. However, Runge-Kutta requires four force calculations per time step, while SPIRAL only requires one. Although the error vs the time step provides insight into the accuracy of each algorithm, it is difficult to evaluate their stability. Therefore, to assess the stability of each algorithm, the plot of the error as a function of time using a time step of $$dt = 10^{-5}$$.

Code that reproduces these benchmarks and other benchmarks is publically available.

Non-Leapfrog Version
To derive a non-leapfrog variant of SPIRAL, we follow the same procedure of the derivation section but with $$a=t $$:

$$q(t + \Delta t) = q(t)\,\exp\left[\frac{\Delta t}{2}\omega(t) + \frac{\Delta t^2}{4}\dot{\omega}(t) + \mathcal{O}(\Delta t^3)\right]\,, $$

which is the non-leapfrog version of SPIRAL. Using the exponential of a quaternion:

$$ q(t+\Delta t) = q(t)\left(\cos{\theta_2} + \frac{\vec{\omega}\left(t\right)}{\left|\vec{\omega}\left(t\right)\right|}\sin{\theta_2} \right)\left(\cos{\theta_3} + \frac{\dot{\vec{\omega}}\left(t\right)}{\left|\dot{\vec{\omega}}\left(t\right)\right|}\sin{\theta_3} \right) $$,

with

$$\theta_2 = \frac{\Delta t}{2}\left|\vec{\omega}\left(t\right)\right| $$,

$$\theta_3 = \frac{\Delta t^2}{4}\left|\dot{\vec{\omega}}\left(t\right)\right| $$.

The non-leapfrog variant does not need a half backward step for initialization of the angular velocity. However, in this version the quaternion should be updated before the angular velocity.