Draft:SPIRAL algorithm

• Comment: This is an article about an algorithm published in 2024, and is a combination of a textbook set of notes on it and an advert. Wikipedia is not the place for this. You need to find 12 or so completely independent reviews or references to it first, which will probably not occur fir some years. Ldm1954 (talk) 12:53, 15 May 2024 (UTC)

---

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: [1]

${\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:

using StaticArrays
using ReferenceFrameRotations

function norm(v::Union{AbstractArray, Quaternion})
    return sqrt(mapreduce(x -> x*x, +, v))
end
function Lab_to_body(v::SVector{3, Float64}, q::Quaternion{Float64})::SVector{3, Float64}
    return vect(inv(q)*v*q)
end
function Body_to_lab(v::SVector{3, Float64}, q::Quaternion{Float64})::SVector{3, Float64}
    return vect(q*v*inv(q))
end

function w_dot(w::SVector{3, Float64}, M::SVector{3, Float64}, II::SVector{3, Float64})::SVector{3, Float64}
    return @inbounds SVector{3, Float64}(
            (M[1] + w[2]*w[3]*(II[2] - II[3]))/II[1],
            (M[2] + w[3]*w[1]*(II[3] - II[1]))/II[2],
            (M[3] + w[1]*w[2]*(II[1] - II[2]))/II[3]
        )
end
function SSPRK3(w::SVector{3, Float64}, M::SVector{3, Float64}, II::SVector{3, Float64}, dt::Float64)::SVector{3, Float64}
    # SSPRK3
    k1::SVector{3, Float64} = dt*w_dot(w                 , M, II)
    k2::SVector{3, Float64} = dt*w_dot(w + k1            , M, II)
    k3::SVector{3, Float64} = dt*w_dot(w + 0.25*(k1 + k2), M, II)
    return w + (k1 + k2 + 4.0*k3)/6.0
end
"""
- q: Current particle orientation.
- w: Current angular velocity (in the lab frame) of the particle.
- M: Torque (in the lab frame) acting on the particle.
- II: Inertia tensor in the principal axis frame. In this frame, the tensor is diagonal. Then, we assume it's a 3D vector. 
- dt: Time step.
"""
function step(q::Quaternion{Float64}, w::SVector{3, Float64}, M::SVector{3, Float64}, II::SVector{3, Float64}, dt::Float64)
    w = Lab_to_body(w, q)
    M = Lab_to_body(M, q)

    # Update velocity w(t + dt/2)
    w = SSPRK3(w, M, II, dt)

    # Update quaternion 
    θ1::Float64 = norm(w)*dt/2.0
    q = q*Quaternion(cos(θ1), sin(θ1)*w/norm(w))

    # Go back to lab frame    
    w = Body_to_lab(w, q)
    return (q, w)
end

Derivation

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

${\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[3]

$e^{\theta {\hat {u}}}=\cos {\theta }+{\vec {u}}\sin {\theta }$,

Comparison of the quaternion relative error for different algorithms and SPIRAL at t=1 s as a function of the time step, as compared to an analytical solution.

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

Comparison of the angular velocity relative error for different algorithms and SPIRAL at t=1 s as a function of the time step, as compared to an analytical solution.

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.

Comparison of the quaternion relative error for different algorithms and SPIRAL at t=1 s as a function of simulation time, as compared to an analytical solution.

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

Comparison of the angular velocity relative error for different algorithms and SPIRAL at t=1 s as a function of simulation time, as compared to an analytical solution.

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: [4]

$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.