This submission is not suitable for Wikipedia. Please read
"What Wikipedia is not" for more information.
If you would like to continue working on the submission, click on the "Edit" tab at the top of the window.
If you have not resolved the issues listed above, your draft will be declined again and potentially deleted.
If you need extra help, please ask us a question at the AfC Help Desk or get live help from experienced editors.
Please do not remove reviewer comments or this notice until the submission is accepted.
Where to get help
If you need help editing or submitting your draft, please ask us a question at the AfC Help Desk or get live help from experienced editors. These venues are only for help with editing and the submission process, not to get reviews.
If you need feedback on your draft, or if the review is taking a lot of time, you can try asking for help on the
talk page of a
relevant WikiProject. Some WikiProjects are more active than others so a speedy reply is not guaranteed.
To improve your odds of a faster review, tag your draft with relevant
WikiProject tags using the button below. This will let reviewers know a new draft has been submitted in their area of interest. For instance, if you wrote about a female astronomer, you would want to add the Biography, Astronomy, and Women scientists tags.
Please note that if the issues are not fixed, the draft will be declined again.
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)
Numerical integration algorithm for integrating rigid body equations of motion.
SPIRAL[1] 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.
where is a pure imaginary
quaternion representing the
angular velocity in the reference frame of the rotating body.
Algorithm
SPIRAL[1] applies to both
leapfrog and non-leapfrog architectures, so it's easily adaptable to different code bases. In this context, the
leapfrog version.
The developers of SPIRAL[1] 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
leapfrogalgorithms, the
angular velocity must be initialized half a step backwards. Note that torque is not recalculated for each
Runge-Kutta stage, meaning that SPIRAL[1] only requires one force calculation per time step.
Bellow, a sample implementation in
Julia of a time step:
usingStaticArraysusingReferenceFrameRotationsfunctionnorm(v::Union{AbstractArray,Quaternion})returnsqrt(mapreduce(x->x*x,+,v))endfunctionLab_to_body(v::SVector{3,Float64},q::Quaternion{Float64})::SVector{3,Float64}returnvect(inv(q)*v*q)endfunctionBody_to_lab(v::SVector{3,Float64},q::Quaternion{Float64})::SVector{3,Float64}returnvect(q*v*inv(q))endfunctionw_dot(w::SVector{3,Float64},M::SVector{3,Float64},II::SVector{3,Float64})::SVector{3,Float64}return@inboundsSVector{3,Float64}((M1+w2*w3*(II2-II3]))/II1],(M2+w3*w1*(II3-II1]))/II2],(M3+w1*w2*(II1-II2]))/II3)endfunctionSSPRK3(w::SVector{3,Float64},M::SVector{3,Float64},II::SVector{3,Float64},dt::Float64)::SVector{3,Float64}# SSPRK3k1::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)returnw+(k1+k2+4.0*k3)/6.0end"""- 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."""functionstep(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.0q=q*Quaternion(cos(θ1),sin(θ1)*w/norm(w))# Go back to lab frame w=Body_to_lab(w,q)return(q,w)end
is a pure imaginary
quaternion representing the
angular velocity in the reference frame of the rotating body. Assuming that it depends only explicitly on time and solve by
separation of variables:
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.
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[1], 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[2],
MercuryDPM[3], etc[1]. 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[4],
Johnson algorithm[5],
Fincham Leapfrog[6],
Omelyan advanced leapfrog[7], algorithm found in
MercuryDPM source code [8] and both SPIRAL[1] 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[1], are second-order. SPIRAL[1] is third-order, and
Runge-Kutta4 is fourth-order. However,
Runge-Kutta requires four force calculations per time step, while SPIRAL[1] 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 .
Code that reproduces these benchmarks and other benchmarks is publically
available.[1]
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 :
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.