Time integration of rigid bodies modelled with three rotation parameters

Three rotation parameters are commonly used in multibody dynamics or in spacecraft attitude determination to represent large spatial rotations. It is well known, however, that the direct time integration of kinematic equations with three rotation parameters is not possible in singular points. In standard formulations based on three rotation parameters, singular points are avoided, for example, by applying reparametrization strategies during the time integration of the kinematic equations. As an alternative, Euler parameters are commonly used to avoid singular points. State-of-the-art approaches use Lie group methods, specifically integrators, to model large rigid body rotations. However, the former methods are based on additional information, e.g. the rotation matrix, which must be computed in each time step. Thus, the latter method is difficult to incorporate into existing codes that are based on three rotation parameters. In this contribution, a novel approach for solving rotational kinematics in terms of three rotation parameters is presented. The proposed approach is illustrated by the example of the rotation vector and the Euler angles. In the proposed approach, Lie group time integration methods are used to compute consistent updates for the rotation vector or the Euler angles in each time step and therefore singular points can be surmounted and the accuracy is higher as compared to the direct time integration of rotation parameters. The proposed update formulas can be easily integrated into existing codes that use either the rotation vector or Euler angles. The advantages of the proposed approach are demonstrated with two numerical examples.


Introduction
The parametrization of spatial rotations is an important issue in many fields such as multibody dynamics (MBD) or in spacecraft attitude determination, for example, since the efficiency of non-linear computations involving large rotations depends largely upon the adequacy of the set of parameters adopted [16]. A commonly used type of rotation parameters in multibody systems dealing with the orientation of rigid or flexible bodies are Euler angles [41], whereas in spacecraft attitude determination, the rotation vector components are often used as rotation parameters [39,52]. It is well known, however, that the direct time integration of kinematic equations with three rotation parameters is not possible in singular points. For a detailed discussion of singular points, see [49]. Hence, the rotational motion cannot be reconstructed by direct integration of the kinematic equations.
This contribution provides a novel numerical solution for the problem of direct time integration of Euler angles and, alternatively, the rotation vector as well as a novel implementation strategy for Lie group time integration methods.

The problem with the direct time integration of rotation parameters
To illustrate the problem of direct time integration of rotation parameters, we consider Euler's equations of motion (EOM) for a single rigid body expressed in a local frame [35,40] Jω + ω × (Jω) = τ (1) that provides the angular velocity ω ∈ R 3 , defined by where R is the rotation matrix transforming the body fixed coordinates of a point to the spatial coordinates. In Eq. (2), ω is the skew-symmetric matrix such that ω × x = ωx for ω, x ∈ R 3 [30]. Solving Eq. (1) for angular velocity does not reveal the rotational motion of the rigid body. Therefore, Euler's EOMs (1) must be complemented with the kinematicreconstruction equationṘ which relates the angular velocity ω expressed in a local frame with the time derivative of the rotation matrix [51]. In order to determine R and thus obtain the orientation of the body, it is required to solve Eq. (3), which is an ordinary differential equation (ODE) on the rotation group SO(3) [30,51]. In many applications, however, e.g. rotor dynamics [38], the direct time integration in the rotation group is avoided by the introduction of a parametrization with generalized coordinates 1 q ∈ R 3 , e.g. Euler angles, such that their time derivativesq are related via ω = G(q)q (4) to the angular velocity ω [30]. When solving the kinematic equations (4) numerically, it is important to note that for each choice of rotation parameters q, regardless of the transformation itself, the transformation (4) depends on the configuration. For any three rotation parameters (e.g. Euler angles, Rodrigues parameter, rotation vector, etc.), G becomes singular at certain configurations, leading to discontinuities in q. If such points are passed, the direct time integration of Eq. (4) is no longer possible. Therefore, a reconstruction of the rotational motion is no longer possible, although the discontinuities in q are finite.

State-of-the-art solutions for direct integration of rotation parameters
In spacecraft attitude determination, switching algorithms are often used as a reparametrization strategy to avoid singular points when using three rotation parameters. For example, the singular points of the rotation vector are often avoided by considering only small rotations (thus the rotation vector represents a locally non-singular three-dimensional rotation parametrization) and by accumulating the incremental attitude in a unit quaternion or direction cosine matrix (DCM) such that the rotation angle corresponding to a rotation vector never gets large [39]. However, the storage in a unit quaternion or a DCM increases the amount of memory required as compared to three rotation parameters. In [48] a switching algorithm is proposed to switch between different Euler angle sets to avoid the singular points while integrating the kinematic equations (4). In [29] a reparametrization strategy based on Shuster's idea of artificial rotation [45] is presented, which shifts the rotation angles associated with singular points of the rotation vector in such a way that either the Rodrigues parameters or the rotation axis and rotation angle can be used without running into a singular point. However, applying a switching algorithm during the integration of Eq. (4) increases the computational complexity as well as the implementation complexity compared to the direct integration of Eq. (4). Three rotation parameters can be used in simulation if a reparametrization strategy such as, for example, the updated Lagrangian point of view [9,16] is implemented to avoid singular points [8]. In the updated Lagrangian point of view, rotations are treated in an updated incremental form. In each step of the solution process, the rotation increment required for moving from the known reference configuration to the current configuration is determined by means of the rotation vector, which is used as a local but not as a global parametrization of the rotation group [9]. Using the rotation vector as well as the rotation matrix from the previous time step, the rotation matrix is updated in the updated Lagrangian point of view in each time step. As in the updated Lagrangian point of view the rotation matrix is used as rotation variable instead of three rotation parameters, it may be difficult to implement the updated Lagrangian point of view approach in existing codes that are based on three rotation parameters.
A common approach in MBD and spacecraft attitude determination to avoid singular points is to use Euler parameters (also known as unit quaternions) [38], as they have only four components -just one more parameter than the dimension of the space of rotations. However, the extra component is bothersome, it requires extra storage as well as extra work to manage the redundancy, e.g. renormalization to unit length [17], or requires solving differential-algebraic equations (DAEs) [5,51]. The need of solving DAEs, however, can be eliminated as shown in [51]. However, it is important to note that if the quaternion integration scheme proposed in [51] is implemented correctly, redundancy will not appear at any time.
The rest of the paper is organized as follows: In Sect. 2, we introduce the standard approach to integrating the kinematic-reconstruction equation (3) using Lie group methods. This section presents the fundamentals on which the proposed approach builds. Subsequently, in Sect. 3, we present our approach independently of a specific choice of rotation parameters. In Sect. 4 (resp. Sect. 5), we apply the proposed approach to the rotation vector (resp. Euler angles). The implementation of the proposed approach is discussed in Sect. 6. In Sect. 7, numerical examples are discussed which illustrate the benefits of the proposed approach. Finally, Sect. 8 draws conclusions from the study.

Standard Lie group approach
The singular points in formulations based on three rotation parameters can be circumvented by using the so-called Lie group integrators [3]. In MBD (see, e.g. the integrators proposed in [7,9]), one seeks, similarly as in the Munthe-Kaas methods [18,25,34] for solving ODEs on Lie groups, a solution to the kinematic-reconstruction equation (3) in the form where ∈ so (3) is the new unknown and the solution oḟ In Eq. (5), R 0 is the initial attitude, and so (3) is the Lie algebra of the Lie group SO (3). Note that 0 in Eq. (6) is an initial condition and that so (3) is isomorphic to R 3 ; see, e.g. [3]. The smooth coordinate map φ takes so(3) to SO(3) [22,24]. The most important coordinate map is the matrix exponential 2 exp(·), which may be evaluated on SO(3) efficiently by Rodrigues' formula [3,24]. In Eq. (6), T −1 φ is the inverse tangent operator associated with the coordinate map 3 φ. In case φ( ) = exp( ), is termed instantaneous or incremental rotation vector in the context of Lie group methods on SO(3) [32,47,51]. To ensure a consistent notation in the following, the term incremental rotation vector is associated with regardless of the choice of a coordinate map. The crucial feature of Eq. (5) is that no global parametrization is required. In contrast to Eq. (4), Eq. (6) can be integrated without running into a singular point. This is because serves as a coordinate vector that locally parametrizes the configuration increment from time step i to time step i + 1 and is assumed to describe small rotations [32].
As shown in [7], R 0 is updated in Eq. (5) in each integration step. Therefore, the rotation matrix R 0 represents a kind of history variable that needs to be available in each integration step during the simulation. 4 Handling the history variable R 0 during time integration is cumbersome and requires additional memory compared to three rotation parameters. The latter procedure may therefore be difficult to implement in existing codes. Specifically on embedded systems, additional history variables might be difficult to be integrated.

Proposed approach (independent of specific rotation parameters)
In this section we present the proposed approach for computing consistent updates for a general set of three rotation parameters, which is subsequently applied to the rotation vector as well as Euler angles. 2 In Appx. A.1 and Appx. A.2, we provide the exponential map as well as a first order Cayley transform on SO (3) and the corresponding inverse tangent operators. 3 As shown in [10,36], other coordinate maps than the exponential map or the Cayley transformation can also be used. 4 Note that only the matrix R 0 belonging to the previous time step must be stored. Fig. 1 Schematic comparison of the procedure for computing three rotation parameters. Note that f (R) is a function that computes three rotation parameters from a rotation matrix R In the proposed approach, the rotation matrix R that is reconstructed with Eq. (5) using the incremental rotation vector is parametrized by three rotation parameters 5 q ∈ R 3 , and therefore Eq. (5) may be represented in the form where q 0 represents initial attitude such that R(q 0 ) = R 0 . The idea of the proposed approach is to consider Eq. (7) as that rotation, which is equivalent to a rotation R(q 0 • q) due to an incremental change of the rotation parameters q ∈ R 3 such that applies and therefore to consider as the solution of Eq. (3) instead of Eq. (5). In Eqs. (8)-(9), the symbol '•' represents the composition rule for the selected rotation parameters q. Note that the explicit form of the composition rule '•' for rotation parameters depends on the selected rotation parameters. Note also that the symbol '•' is used here to indicate that the rotation parameters q are determined from the rotation parameters q 0 and q using the appropriate 6 composition rule. In Eq. (9), we assume that and that q does not need to be small. Furthermore, we consider the incremental change of the rotation parameters q as a finite increment of rotation. To be able to compute the rotation parameter update with Eq. (9), q(q 0 , ) must be determined for the selected rotation parameters depending on the incremental rotation vector , where the incremental rotation vector is determined in each time step by solving Eq. (6) numerically. The rotation parameters q can, of course, also be extracted from the rotation matrix as shown in [19,27,43], which is a detour compared to the method proposed in this paper. In the proposed approach, the rotation parameters q are computed directly using the incremental rotation vector without having to evaluate in advance the rotation matrix, see Fig. 1, last line.
As shown in Fig. 1, in the proposed approach the rotation parameters are not determined by integrating the corresponding kinematic equations (4), but by using the incremental rotation vector. Thus, we do not need to deal with the singular points occurring in the kinematic equations (4). Of course, there is no three-parametric description of spatial rotations that is free of singular points [49]. Consequently, singular points are also present in the proposed approach. However, Eq. (9) allows dealing with singular points at the coordinate level. At the coordinate level, singular points can be embedded within the composition rule of rotations that applies to the selected rotation parameters using case distinctions, for example, such that the time integration scheme reconstructs the rotation parameters associated with a singular point without failing.
We would like to note, however, that the idea of computing rotation updates based on the incremental rotation vector is not new at all, see the references in Sect. 1.2 and Sect. 2. However, the authors are not aware that an approach similar to the presented one has already been discussed for the time integration of three rotation parameters.
As a final remark, we would like to note that the reason why we consider different parametrizations for the map φ in Eq. (7) is that the choice of the parametrization of the map φ has a significant influence on the computational effort required to compute q and thus on the computational effort required to determine a rotation parameter update.
The advantages of the proposed approach are shown in the following by the example of the rotation vector, as well as by the example of Euler angles.

Proposed approach for the time integration of rotation vector
In this section, we apply the general approach to integrating rotation parameters presented in Sect. 3 to the rotation vector parametrization, i.e. q = v.

Definition of the rotation vector
If n ∈ R 3 is a unit vector along the rotation axis and ϕ ∈ R is the rotation angle, then the corresponding rotation vector v ∈ R 3 is given by [4,28,32,43] (see also Fig. 2 where the rotation angle ϕ and rotation axis n can be recovered from rotation vector v by

Rotation vector update
To be able to compute an update of the rotation vector in the form of Eq. (9), we need to find the composition rule for rotation vectors depending on the incremental rotation vector 7 .

Fig. 2
Body-fixed vector in positions r before and r * after the rotation (ϕ, n), see also [53] The relation between the incremental rotation vector and the rotation vector The composition rule of two successive rotations with rotation vectors v = ϕ n (14) and with v being the resultant rotation vector, is determined by using Eq. (11) as In Eq. (16), ϕ(ϕ 0 , ϕ, n 0 , n) and n(ϕ 0 , ϕ, n 0 , n) are given by after elementary transformations of the equations given in [1,37,39,54] where I is a 3×3identity matrix. In Eqs. (17)-(18), we have used the relationships and In Eqs. (17)-(18), ϕ 0 • ϕ represents the composition rule for rotation angles and n 0 • n represents the composition rule for rotation axes. It should be noted, however, that Eqs. (17)- (18) can only be evaluated for 8 Equation (13) can be verified with To find a relation between v and , a Cayley transform could also be used as a coordinate map in Eq. (13) instead of the exponential map. However, the relationship between v and would be more complicated since R(v) = cay(v).

Improved rotation vector update
The case = 0 can be covered using the cardinal sine function which is continuous and computable at x = 0 [32]. Therefore, we replace Eqs. (17)- (18) with To cover the case ϕ ∈ S v , we introduce the rotation axisn using the case distinction 9 By using Eqs. (23) and (25), the composition rule for rotation vectors reads which can be evaluated for all ϕ ∈ R without restrictions. The improvements of Eqs. (23) and (24) (resp. Eq. (25)) compared to Eqs. (17)-(18) are shown with the following considerations: • Case 1. In the case that ϕ 0 = 0 and = 0, after a few steps of calculation, it follows from Eqs. (23) and (25) that andn For this case, Eq. (26) therefore gives • Case 2. In the case that ϕ 0 = 0 and = 0, after a few steps of calculation it follows from Eqs. (23) and (25) that andn For this case, Eq. (26) therefore gives • Case 3. Of course, for the trivial case ϕ 0 = 0 and = 0, Eq. (26) yields Note that the rotation vector exhibits a redundancy at the rotation angles ϕ ∈ S v [20]. However, the Rodrigues' formula resolves the redundancy ϕ ∈ S v , Note that the redundancy of the rotation vector at ϕ ∈ S v can be bypassed by scaling the rotation vector such that its magnitude always belongs to an interval ϕ ∈ [−π, π] as shown in [20].
Since we use the exponential map Eq. (77) in Eq. (13) to establish a relation between and v, the inverse tangent operator corresponding to the exponential map T −1 exp in Eq. (82) needs to used in Eq. (6). Therefore, Eq. (6) takes the forṁ The incremental rotation vector is then determined in each time step by solving Eq. (36) numerically.
As a final remark, we would like to note that Engø [14] presented a closed form of the Baker-Campbell-Hausdorff formula on so(3) and has interpreted it as a composition of two rotations, which was extended in [12] for SE(3) and dual quaternions. However, the formula presented in [14] cannot be applied directly for the time integration of the rotation vector, since it only provides correct answers for rotation angles smaller than π /2 as noted in the paper.

Proposed approach for the time integration of Euler angles
In this section, we apply the general approach to the Euler angle parametrization, i.e. q = α.
To be able to compute an update for Euler angles in the form of Eq. (9), we need to find the composition rule as well as the Euler angle increment depending on the incremental rotation vector . Since there are 12 possible sets of Euler angles [13,19,44], we present a general procedure to determine the Euler angle increment for an arbitrary set of Euler angles. In Appx. B we apply the general procedure shown in this section to determine the Euler angle increment for Cardan/Tait-Bryan angles, which has been omitted here for reasons of space and clarity.
Any combination of three elementary rotations is given as with e i being a constant unit vector [32,48] and exp being the exponential map Eq. (77). For the Cardan/Tait-Bryan angle parametrization, these are the unit vectors whereas for 313-Euler angles these are unit vectors [32] In the proposed approach, the composition rule of two successive rotations with Euler angles α 0,1 , α 0,2 , α 0,3 and α 1 , α 2 , α 3 where α 1 , α 2 , α 3 are the resulting Euler angles is applied as 10 where α i is implicitly defined by the standard Lie group method Eq. (5) and explicitly computed in our approach. By using Eq. (40), the elementary rotation exp(α i e i ) associated with the ith Euler angle α i is written as By multiplying Eq. (41) with exp T (α 0,i e i ) from the left, the elementary rotation exp( α i e i ) associated with the Euler angle increment α i is given by From Eq. (42) two equations are obtained for α i , one equation for the sine and one equation for the cosine of an incremental angle, which read 11 cos α i = cos α 0,i cos α i + sin α 0,i sin α i =: From Eqs. (43)-(44), the Euler angle increment α i is computed using the proper fourquadrant inverse tangent function atan2 12 In order to be able to compute α i with Eq. (45), sin α i and cos α i in Eqs. (43)-(44) need to be expressed depending on the incremental rotation vector. However, the equations that provide this connection assume different forms for each of the 12 possible Euler angle 10 Note that the Euler angle composition rule shown in Eq. (40) is written for a single elementary rotation whose axis of rotation is given by one of the unit vectors from Eq. (38) (resp. Eq. (39)). Note also that α i is associated with the ith Euler angle increment and that α i is non-linear in α 0,1 , α 0,2 , α 0,3 and , see Eq. (45). 11 Equations (43)-(44) always take the same form, independent of the particular Euler angle set, which can be verified by direct calculation of Eq. (42) with different combinations of e i . 12 In the numerical examples shown in Sect. 7, we have used the atan2 function implemented by default in Matlab2018b. We would like to point out that compared to the implementation of the atan2 function proposed in IEEE [21], the Matlab2018b implementation of the atan2 differs by having the property atan2(0,0)=0. This property is essential for the correct implementation of the proposed Euler angle integration scheme.
combinations. The procedure for establishing a relationship between sin α i and cos α i and is illustrated in Appx. B using the Cardan/Tait-Bryan angles as an example. By combining the Euler angles α 0,i in vector and the Euler angle increments α i obtained with Eq. (45) in vector and the Euler angles α i in vector the composition rule of two successive rotations in vector form for Euler angles is

Time integration scheme for the proposed approach
In the proposed approach, the rotational update within the ith integration step is based on Eq. (9), where the incremental rotation vector updates a rotation of q i to q i+1 . Therefore we express the update for the step i + 1 in the form where i is the incremental rotation vector of step i and q i represents the rotation parameters at the ith step; see, e.g. [3,8]. In order to obtain i , we write Eq. (6) aṡ where 0 is the initial condition for the integration of Eq. (51) in the ith integration step and ω i+1 is the angular velocity of step i + 1. The equations (51) should be integrated within each integration step together with the dynamical equations of motion, which determine the angular velocity field ω from the angular acceleration fielḋ by using Euler's EOMs (1); see, for example, [51].

Explicit Euler (RK1)
Using the time step size h, the explicit Euler method takes the following form for the proposed approach: which differs from the standard explicit Euler Lie group method presented in [26] by the fact the no history variables are required. The conventional explicit Euler method for the direct time integration of rotation parameters would have the form The comparison of Eqs. (56)-(57) with Eqs. (53)-(55) shows that the proposed approach can be integrated into existing applications that use three rotation parameters by simply replacing 13 Eq. (57) with Eqs. (54)-(55). Using the rotation vector, i.e. q = v, for Euler's EOMs (1) for the proposed approach, the explicit Euler method takes the form At this point, we would like to give a brief interpretation of the RK1 method with regard to its functionality in a singular point: Sinceω is bounded (due to the physical limitation of torque), ω is also bounded, and thus also . Consequently, v is also bounded due to the definition of the composition rule for rotation vectors defined in Sect. 4. However, if one would consider the direct integration of the kinematic equations of the rotation vector Eq. (87) (Appx. A.3),v would not be bounded at a singular point.
Note that, according to the explicit Euler method, Eqs. (53)-(55) are only first-order accurate and have a numerical stability limit for the step size h in case of stiff differential equations.

Explicit Runge-Kutta method of fourth-order (RK4)
To achieve a higher accuracy with larger step sizes h, we consider the explicit fourth-order Runge-Kutta method presented for time integration of Euler parameters by Terze et al. in [51] in all presented numerical examples, which we have written here in accordance with the proposed approach. Starting with values at the previous step q i , ω i , the slope estimations are obtained by The four coefficients K 1 , . . . , K 4 represent the slopes within one time step (i → i + 1) in the solution process of Eq. (51) and the four coefficients k 1 , . . . , k 4 contain the slopes within one time step in the solution process of the dynamic equations of motion. The equations update the angular velocity vector ω and the rotation parameter vector q to (q i+1 , ω i+1 ) based on the current state (q i , ω i ).
Since the composition rule of two successive rotations can be expressed in a closed form for the rotation vector and the Euler angles, the order of accuracy of the overall algorithm depends only on the accuracy of the ODE integrator used to solve Eq. (51). The main paradigm of the Lie group methods is the reformulation of the underlying equation on the Lie group as an algebra action, and since the Lie algebra is a linear space, all reasonable discretization methods can be expected to respect its structure [23,24]. Therefore the user can influence the accuracy of the overall algorithm by choosing the ODE integrator to solve Eq. (51).

Behaviour in singular configurations
As a first numerical example, we investigate the behaviour of the proposed rotation vector and Euler angle approach near and in a singular point where G(q) in Eq. (4) becomes singular. For this purpose, we consider a parameter study for the torque-free, free rotation of a rigid body with parametrized initial conditions, which are changed such that the orientation approaches the singular point and finally reaches it. Since there are a total of 12 possible Euler angle sets, we limit ourselves to Cardan/Tait-Bryan angles for the proposed Euler angle approach in all following numerical examples. The rigid body is considered to be a box with inertia tensor J given in standard units equal to J = diag (5.2988, 1.1775, 4.3568). The rotation vector parametrization has a singular points at the rotation angles ϕ ∈ S v , see Eq. (87) in Appx. A.3, whereas the Cardan/Tait-Bryan angles exhibit a singular point at α 2 = ±n π 2 with n ∈ N, see Eq. (88) in Appx. A.4. In order to determine ω(t), (t) and thus also v(t) and α(t), Euler's equations of motion (1) as well as the incremental rotation vector differential equation (36) are integrated using the RK4 algorithm presented in Sect. 6.2, whereby in the proposed Cardan/Tait-Bryan approach we use the exponential map Eq. (78) as a coordinate map. Similar results for Cardan/Tait-Bryan angles can be obtained by choosing the Cayley transform Eq. (83) as a coordinate map. The whole algorithm has been implemented in MATLAB R2018b.

Numerical tests for the proposed rotation vector approach
In this subsection we investigate the behaviour of the rotation vector approach proposed in this paper near and in the singular point at ϕ = 0. Similar results can be obtained for other singular points of the rotation vector. In order to reach the singular point at ϕ = 0 step by step, we change the initial condition for angular velocity according to gradually increasing the factor ∈ {0, 1e-7, 1e-5, 1} whereby the initial orientation of the rigid body is kept constant at v 0 = 0, − π 2 , 0 T rad. Thus, the singular point is reached after t = 0.5 s. As can be seen in Fig. 3(a), in the case of ϕ = 0, the rotation vector v StdAppr can no longer be determined by integrating Eq. (87) at t = 0.5 s due to the division by zero. In contrast, as can be seen in Fig. 3(b), the proposed approach for computing a rotation vector update v PropAppr is able to pass through the singular point without failing. Figure 4 illustrates convergence in the norm of the rotation error v Ref − v converged for decreasing values of the integration step h n = (2 (1−n) ) n=2,3,..., 15 , where the rotation vector v StdAppr is obtained by directly integrating (87). The converged reference solution v Ref is obtained with h = 1/32768 s using the proposed approach. As can be seen in Fig. 4(a), the rotation vector v StdAppr converges slower the closer the rotation angle comes to the singular point at ϕ = 0 and it even fails as the rotation angle reaches ϕ = 0. The failure of the standard method corresponds to the expectation, since in the kinematic Eq. (87) a division Fig. 4 Error in the norm of the rotating rigid body's rotation vector shown on a double-logarithmic scale by zero occurs at ϕ = 0. In contrast, the rotation vector v PropAppr computed with the proposed approach exhibits the expected fourth-order convergence characteristics for all considered values of , as can be seen in Fig. 4(b). Furthermore, as can be seen in Fig. 4

Numerical tests for the proposed Cardan/Tait-Bryan angle approach
In this subsection we investigate the behaviour of the proposed Cardan/Tait-Bryan angle approach near and in the singular point at α 2 = π 2 . In order to reach the singular point at α 2 = π 2 step by step, we change the initial condition for angular velocity according to gradually increasing the factor ∈ {0, 1e-5, 1e-2, 1e-1} whereby the initial orientation of the rigid body is kept constant at α 0 = 0, 0, 0 T rad. Thus, the singular point is reached after t = 0.5 s. As can be seen in Fig. 5(a), as α StdAppr 2 (which is dominated by the rotation angle α 2 ) comes to the singular point at α 2 = π 2 , α StdAppr 2 can no longer be determined at t = 0.5 s by integrating the kinematic equation (88). This corresponds to the expectation, since in Eq. (88) a division by zero occurs at α 2 = π 2 . In contrast, as can be seen in Fig. 5(b), the proposed approach is able to pass through the singular point at α 2 = π 2 without failing. In the case that = 0, the rigid body performs a pure rotation around the middle axis of rotation. For this simple case an analytical solution of Eq. (88) for the rotation angle α 2 can be given, which is However, as can be seen in Fig. 6(b), the Cardan/Tait-Bryan angles determined with the proposed approach show finite jumps compared to the analytical solution ( Fig. 6(a)) at angle values for α 2 = n π 2 with n ∈ N (these are the angle values where the Cardan/Tait-Bryan angle parametrization has singular points). However, the Cardan/Tait-Bryan angles determined with the proposed approach can, if desired, be converted in the course of post-processing, see Fig. 7(a). The conversion can be done by evaluating the rotation matrix and extracting the angles as shown for example in [19]. However, we would like to emphasize that this conversion is not necessary in order to uniquely represent the rotational attitude of a rigid body using the proposed approach. Figure 7(b) shows the deviation of the converted angle α 2 from the analytic solution for a time interval of 100 s, whereby the angles determined with the proposed approach were determined with a time step size of h = 1e-3 s.   Fig. 8(a), α StdAppr converges slower the closer α 2 comes to the singular point at α 2 = π 2 and it even fails as the rotation angle α 2 reaches α 2 = π 2 . The failure of the standard method corresponds to the expectation, since in the kinematic Eq. (88) a division by zero occurs at α 2 = π 2 . In contrast, the proposed Cardan/Tait-Bryan angle approach exhibits the expected fourth-order convergence characteristics for all considered values of , as can be seen in Fig. 8(b). Furthermore, as can be seen in Fig. 8(b), α PropAppr 2 is in case that α 2 = π 2 already converged at larger time steps. For the convergence investigations in Fig. 8, the Cardan/Tait-Bryan angles determined with the proposed approach were converted as described above in the case of = 0.

Influence of the used ODE integrator for solving Eq. (51)
To illustrate how the accuracy of the proposed approach depends on the order of the ODE integrator used to integrate Eq. (51), we compare in Fig. 9 the rotation vector error when the integration of Eq. (51) (resp. Eq. (36)) was performed using the two methods RK1 (see Sect. 6.1) and RK4 (see Sect. 6.2). As model problem we consider here the rotation of a torque-free rigid body about an axis close to its unstable axis of rotation. The rigid body is considered to be a box with inertia tensor equal to J = diag(5.2988, 1.1775, 4.3568). Since  Convergence is investigated at t = 0.5 s. As can be seen in Fig. 9, the overall integration results show first and fourth order convergence characteristics as expected. This demonstrates the fact that the order of accuracy of the proposed approach can be selected by the user and that it depends only on the accuracy of the ODE integrator used for the integration of Eq. (51). The versatility of the proposed approach in terms of the possibility of arbitrarily selecting the order of accuracy according to the field of application should contribute positively to its use in the various fields of application, including those requiring higher-order integration schemes. The same behaviour is observed in the Euler angle integration scheme presented in this paper. However, due to space reasons we do not present the corresponding figure here. We would like to note that the latter behaviour is also observed for the Euler parameter integration scheme presented in [51].

Convergence comparison with the standard Lie group method
To illustrate how the proposed approach compares to the standard Lie group method in terms of accuracy, we compare in Fig. 10 Table 1 shows the position error p Ref − p converged of the three methods considered.
As can be seen in Fig. 10 and Table 1, the proposed approach exhibits almost the same convergence behaviour as the standard Lie group method. The difference between the standard Lie group method and the proposed approach that can be observed in Table 1 is due to round-off errors. From this observation we conclude that the proposed method is equivalent to a Lie group method. As the proposed approach is equivalent to a Lie group method, we will not further compare the proposed approach with the standard Lie group method in the following accuracy comparisons. For further investigations and results on Lie group methods, we would like to refer to the corresponding literature; see, e.g. [2,3,18,26].

Heavy top
As a second example, we present the integration of the dynamics of a heavy top to illustrate the performance of the proposed approach. The heavy top is a kind of benchmark problem in the context of geometric integration methods, which is, for example, also discussed in [3,8,15,33,50,51]. Following [51], the configuration space of the heavy top is set as SO (3) and the dynamical model is formulated in the classical ODE form on the basis of Euler's rotational equation In Eq. (74), ω andω are the body angular velocity and angular acceleration, J FP is the tensor of inertia with respect to the top's fixed point, which is computed from the mass center inertia tensor J using the parallel-axis theorem In Eq. (75), m is body mass and r b is the body mass center position expressed in the body fixed frame with respect to the top's fixed point. The gravity torque τ FP with respect to the fixed point is computed as where g = [0, 0, −9.81] T represents gravity and R is a rotation matrix. 14 In standard units, the mass center inertia tensor is defined as In order to determine ω(t), (t) and thus also the rotation vector v(t) as well as the Cardan/Tait-Bryan angles α(t), the equations of motion (74) as well as the incremental rotation vector ODE (36) are integrated using the RK4 algorithm shown in Sect. 6.2. In the proposed Cardan/Tait-Bryan angle approach, we use the exponential map as a coordinate  Figs. 11-13 and 17 can, of course, also be achieved with the proposed Cardan/Tait-Bryan angle approach, which is why we do not present them explicitly here.
The time evolution of the rotation vector components as well as the time evolution of the Cardan/Tait-Bryan angle components, obtained by the proposed approach, are illustrated in Figs. 15-16. The deviation of the orthogonality constraint of the rotation matrix is shown in Fig. 17(a). As can be seen in Fig. 17(a), the proposed approach does not introduce any drift in the integration results during the whole simulation time domain of 1000 s. In comparison, Fig. 17(a) shows the deviation of the orthogonality condition of the rotation matrix in the case where Euler parameters are used as rotation parameters, whereby the Euler parameters were com-  [51]. As illustrated in Fig. 17(b), the maximum error of the orthogonality constraint of the rotation matrix is for the simulation time domain of 1000 s in the range of machine precision. Figure 18 illustrates the convergence in norm of the rotation error v − v converged as well as the convergence in norm of the rotation error α − α converged 2 for decreasing values of the time step h n = (10 −2 · 2 (1−n) ) n=1,2,..., 11 . The reference solution for the rotation vector v converged = v(t = 1) as well as the reference solution for the Cardan/Tait-Bryan angles α converged = α(t = 1) were computed using the proposed approach and a integration step of h = 1/204800. In order to be able to compare the convergence behaviour of the proposed approach with the convergence behaviour of the standard approach (direct integration of the kinematic equations of the rotation vector and the Cardan/Tait-Bryan angles), the initial orientation of the heavy top may not be defined by v 0 = 0, since the kinematic equations (87) of the rotation vector cannot be evaluated at v 0 = 0 (because of the divi- sion by zero). Therefore, the initial conditions for the convergence considerations are set as As can be seen in Fig. 18(a) as well as in Fig. 18(b), the proposed approach for computing a rotation vector update and a Cardan/Tait-Bryan angle update exhibit smaller errors (of about the order of 1e-2) in norm of the rotation error for each time step h in comparison to the rotation updates computed by directly integrating the corresponding kinematic equations. Especially the newly introduced rotation vector formulation shows a much smaller error at larger time steps than the direct integration of the corresponding kinematic equations, as can be seen in Fig. 18(a).

Computational costs
Finally, we compare in Tables 2-3 the number of arithmetic operations required by the proposed approach Eq. (9) with those required by the standard Lie group method Eq. (5) and Fig. 17 Time evolution of the norm of rotation matrix orthogonality condition as well as maximum error in rotation matrix orthogonality condition computed for one million steps by the Euler parameter integration scheme presented in [51] to determine an update for the rotation vector or Euler angles. For the comparison of the standard Lie group method and the Euler parameter integration scheme with the proposed approach, three rotation parameters are computed from the rotation matrix Eq. (5), respectively the Euler parameters computed with Eq. (12) in [51]. The rotation matrix R 0 in Eq. (5) and the Euler parameters θ 0 are reconstructed from the rotation vector v 0 and the Cardan/Tait-Bryan angles α 0 , respectively. The cost of the integration of the incremental rotation vector differential equation (6) is left out in the following considerations, since the computational cost depends on the chosen integration scheme. As the three rotation parameters determined by direct integration of the kinematic equations (4) cannot be determined in singular points, we do not consider the computational costs of this methods any further here. The columns in Tables 2-3 show the respective number of additions, subtractions, multiplications, divisions, trigonometric functions and roots and the results of adding all these together. While the rows of Tables 2-3 show the method considered in each case.
As shown in Table 2, the proposed rotation vector approach requires fewer operations to determine a rotation vector update as compared to determining the rotation vector update from the rotation matrix computed with the standard Lie group method. In contrast, the method which computes an update for the rotation vector from the Euler parameters requires fewer operations than the proposed approach. As shown in Table 3, the proposed  Table 2 Number of arithmetic operations needed in the proposed approach Eq. (26), in the standard Lie group method Eq. (5), and in the quaternion integration scheme [51] for computing a rotation vector update. Note that f (·) is a function that computes the rotation vector from a given rotation matrix R according to page 54 in [46], respectively Euler parameters according to page 482 in [42]. Note also that the term exp or cay in the second and third row indicates which coordinate map is used in Eq. ( Cardan/Tait-Bryan angle approach requires fewer operations to determine a Cardan/Tait-Bryan angle update as compared to determining the Cardan/Tait-Bryan angle update from the rotation matrix. In contrast, the method which computes an update for the Cardan/Tait-Bryan angles from the Euler parameters requires fewer operations as the proposed approach.

Conclusions
A novel numerical solution for solving rotational kinematics by using three rotation parameters is presented in this paper. The proposed approach is illustrated by the example of the rotation vector and the Euler angles. In contrast to standard formulations based on three rotation parameters, in which singular points are avoided, for example, by applying reparametrization strategies, singular points can be passed through in the proposed approach by computation of finite increments of rotation. On the contrary to standard formulations based on three rotation parameters, the proposed approach is based on the numerical integration of the kinematic relations in the form of the so-called incremental rotation vector. The kinematic relations of the incremental rotation vector form a system of ODEs on the Lie algebra so(3) of the rotation group SO(3) that can be solved singularity-free by any standard ODE integration scheme. In the proposed approach, after the incremental rotation vector has been determined for the current step, the rotation parameter update for the current step is determined with the composition rule of rotations that applies to the particular rotation parameters using the finite increments of rotation (which are expressed in dependence of the incremental rotation vector), and the rotation parameters of the previous integration step. By taking this route, the proposed approach is able to reconstruct the rotation parameters even in singular configurations. As shown in the paper, the proposed approach exhibits the same accuracy as the standard Lie group method and the accuracy is higher as compared to the direct time integration of rotation parameters. As the presented results show, the order of accuracy of the presented integration scheme depends only on the accuracy of the ODE integrator used for solving Eq. (51). Therefore the user can influence the accuracy of the overall algorithm by choosing the ODE integrator to solve Eq. (51).
The new methods proposed in this paper are more favourable compared to the standard Lie group method in terms of their computational effort. However, the newly proposed methods are less favourable in terms of the required computational effort as compared to reconstructing the rotation vector or Cardan/Tait-Bryan angles from the Euler parameter integration scheme proposed in [51]. Therefore, the improvement of the proposed approach regarding its computational efficiency will be the subject of future research.
For solving stiff problems, the use of an implicit time integration scheme is advantageous [18]. Therefore, embedding the proposed approach in an implicit time integration scheme will be the content of future work.
In summary, the current paper provides a novel numerical solution for the long-standing problem of direct time integration of Euler angles and, alternatively, the rotation vector. The paper provides ready to implement update formulas that allow to compute consistent updates for the rotation vector and the Euler angles. The provided update formulas can be easily integrated into existing codes that uses either the rotation vector or Euler angles and provide a novel implementation strategy of Lie group methods.
The inverse tangent operator T −1 exp corresponding to the exponential map on SO(3) is given by [51] T −1 exp ( ) = I + which exhibits singular points at = ±2πk, k = 0, 1, 2, . . . Decisive for the evaluation of Eq. (80) is the situation = 0, because T −1 exp has a discontinuity there, which is removable due to lim →0 det(T −1 exp ( )) = 1, as shown in [32]. In the implementation, the singular point at = 0 in Eq. (80) is according to Eq. (81) removed by using the case distinction which allows a singular-point-free integration of the angular kinematics in the range ∈ (−2π, 2π), i.e. for a full double turn [32].

A.2 The Cayley transform on SO(3)
The Cayley transform can be used as coordinate map φ( ) = cay( ) as well. This is due the fact that the Cayley transform also maps algebra elements to group elements, which is well know on SO(3) [31].
The Cayley transform cay : R 3 → SO (3) is defined by [10,26] cay( ) = I − while the corresponding inverse tangent operator is given by
From Eqs. (98)-(103), the Cardan/Tait-Bryan angle increments α i can be computed without restrictions using the atan2 function Note that it is not necessary to select a specific coordinate map φ( ) in Eq. (82) when deriving the general expression for the angle increments depending on the incremental rotation vector. Therefore, it is possible to implement the proposed approach for any choice of Euler angles in such a way that it is easy to switch between different coordinate maps. This allows influencing the computational efficiency of the proposed approach in a specific way. Note also that, in general, the coordinate map and the inverse tangent operator must match, i.e. T −1 φ = T −1 exp if φ = exp or T −1 φ = T −1 cay if φ = cay. However, note that switching between different coordinate maps is not necessary in the proposed approach. As a final remark, we would like to note that for an efficient implementation of the proposed Cardan/Tait-Bryan angle approach, it is advantageous to optimize the implementation for a specific coordinate map by computing and implementing u i and v i explicitly.