An improved convergence based on accelerated modulus-based Gauss–Seidel method for interactive rigid body simulations

In this paper, we propose a novel method that fits linear complementarity problems arising in interactive rigid-body simulations, based on the accelerated modulus-based Gauss–Seidel (AMGS) method. We give a new sufficient condition for the convergence of the generated sequence under a milder condition on the matrix splitting than the special case of the AMGS method. This gives a flexibility in the choice of the matrix splitting, and an appropriate matrix splitting can lead to a better convergence rate in practice. Numerical experiments show that the proposed method is more efficient than the simple application of the AMGS method, and that the accuracy in each step of the proposed method is superior to that of the projected Gauss–Seidel method.


Introduction
In rigid-body simulations, interactions (for example, normal forces) between rigid bodies are often mathematically modeled as constraints. For computing these constraint forces, we usually need to solve certain equations. Two representative constraint formulations for rigid-body simulations are acceleration-based formulations [3] and velocity-based formulations [1]. In the acceleration-based formulations, the constraints are described with forces and accelerations of rigid bodies; we first compute forces and accelerations, then integrate them to obtain velocity changes. On the other hand, in the velocity-based formulations, the variables in the constraints are impulses and velocities of the rigid bodies. In this paper, we focus on the velocity-based formulations, since the velocity-based formulations are widely used and are known to be superior to the acceleration-based formulations in many aspects (for example, see [9,14]). In recent years, position-based formulations [8] have also been developed for rigid-body simulations. The position-based method was originally proposed for cloth simulations [7], but it is widely used for various simulations especially in the context of computer graphics [6].
There are two main categories for solving constraints, iterative approaches and direct approaches. Our interest in this paper is the iterative approaches rather than the direct approaches, since the direct approaches such as pivoting methods often suffer from time complexity and numerical instability as pointed in [13]. In the impulse-based iterative approaches [9], impulses are applied to the rigid bodies sequentially, until certain convergence conditions are satisfied. Tang et al. [16] proposed an impulse-based energy tracking method that computes accurate velocities by applying impulses iteratively. This method implicitly computes relative velocities after collisions, since an explicit computation of the velocities may lead to physically inaccurate results when multiple contacts occur simultaneously.
Contact constraints are frequently modeled in a form of complementarity problems [4]; in particular, linear complementarity problems (LCPs) give mathematical formulations for frictionless contacts. For frictional contacts, we can use nonlinear complementarity problems (NCPs) to formulate accurate Coulomb friction [18], or, as we will discuss later in Sect. 4.1, we can use boxed LCPs for approximated contact friction which employs friction pyramids instead of accurate friction cones. Among various iterative methods for solving LCPs, the projected Gauss-Seidel (PGS) method [5] has many extensions for solving contact constraints [9,10,14]. Another iterative approach for solving LCPs is the use of modulus-based methods. Bai [2] established modulus-based matrix splitting iteration (MMSI) methods, which include modulusbased Jacobi (MJ), modulus-based Gauss-Seidel (MGS), modulus-based successive over relaxation (MSOR), and modulus-based accelerated overrelaxation (MAOR) iteration methods as special cases. Mezzadri and Galligani [12] proposed an extension of the MMSI methods so that they can be used to solve horizontal linear complementarity problems (HLCPs). Zheng and Vong [19] examined the convergence of the MMSI methods for HLCPs, and they proposed a more general convergence result. Zheng and Yin [20] proposed accelerated modulus-based matrix splitting iteration (AMMSI) methods as an improvement of Bai [2]. In a similar way to the MMSI methods, the AMMSI methods include accelerated modulus-based Jacobi (AMJ), accelerated modulus-based SOR (AMSOR), and accelerated modulus-based accelerated overrelaxation (AMAOR) iteration methods. Furthemore, the AMMSI methods also devised the accelerated modulus-based Gauss-Seidel (AMGS) method.
In this paper, we give a theoretical proof on the convergence of the accelerated modulus-based Gauss-Seidel (AMGS) method. Zheng and Yin [20] already discussed the convergence in a general case, but their assumption for the case of positive-definite coefficient matrices is too restrictive to apply the same discussion to rigid-body simulations. We propose another sufficient condition of the convergence of the AMGS method when the coefficient matrix of the LCP is positive definite, so that we can choose a parameter which leads to a faster convergence. We also show a simple example that is covered by the condition we propose, but is not covered by the condition of a general case proposed in [20].
We also improve the efficency of the proposed method. In many applications of real-time simulations, interactive computer graphics and operations are considered most important, since, if the computation in each simulation step is considerably slower than real time, the quality of the users' experience would be seriously degraded. Since the AMGS method proposed in [20] is not designed for interactive rigid-body simulations, a simple application of the AMGS method causes inefficiency and it is a serious disadvantage for real-time simulations. The proposed method focuses the update formula in the AMGS method and exploits the structures related to the generalized velocity vector of rigid bodies.
Through numerical experiments, we observed that the proposed AMGS method attained shorter computation time than the original AMGS method. Furthermore, its convergence rate in each iteration was better than that of the PGS method. These results indicate that the proposed method is useful for practical real-time simulations. Mezzadri [11] showed that under a specific condition, the PGS method and the AMGS method are equivalent in that AMGS iterations can be written as PGS iterations. Mezzadri also pointed out that the AMGS method also performs like the projected successive over-relaxation (PSOR) method under a specific parameter choice, and this is consistent with our numerical results.
The outline of this paper is as follows. In Sect. 2, we briefly introduce a formulation of velocity-based constraints as an LCP. We also discuss the PGS method and the AMGS method to solve LCPs. We prove convergence theorems of the AMGS method in Sect. 3, and the application of the AMGS method to rigid-body simulations is developed in Sect. 4. In Sect. 5, we will show numerical results to verify the efficiency of the AMGS method. Finally, we will give a conclusion in Sect. 6.

Linear complementarity problem with velocity-based constraints
For the latter discussions, we briefly introduce an LCP that arises from velocity-based constraints. For more details, the readers can refer to [4,15,17]. During a rigid-body simulation, we keep tracking movements of the rigid bodies in multiple time periods, therefore, an entire simulation is divided into a sequence of simulation steps, and each simulation step corresponds to a small time step. Since the constraints on the rigid bodies should be satisfied at each time, we solve the following LCP in each simulation step: In this paper, we use the superscript T to denote the transpose of a vector or a matrix.
The decision variable in this LCP is ∈ ℝ m , which is the impulse vector applied to the rigid bodies in the constraint space. The first constraint in (1) requires to be nonnegative to ensure that the constraint impulse must be repulsive.
The second constraint in (1) corresponds to the velocity constraints in the rigid-body simulation. The vectors b ∈ ℝ m and v ∈ ℝ n are the bias vector in a constraint space and the generalized velocity vector of rigid bodies, respectively. More precisely, when we have N rigid bodies, v is a vector that consists of N linear and angular velocities, i.e., where v i ∈ ℝ 3 and i ∈ ℝ 3 are the linear velocity and the angular velocity of the ith rigid body, respectively, for i = 1, … , N ; thus the length of v is n = 6N . The matrix J ∈ ℝ m×n is the Jacobian matrix corresponding to the velocity constraints. The generalized mass matrix of the rigid bodies M ∈ ℝ n×n consists of masses and inertia tensor matrices in the diagonal positions: where m i ∈ ℝ and I i ∈ ℝ 3×3 are the mass and the inertia tensor matrix of the ith rigid body, respectively. We also use E r ∈ ℝ r×r to denote the identity matrix of order r. The inertia tensor matricies I 1 , … , I N are symmetric, so is M.
The third constraint in (1) is a complementarity condition. We can understand this complementarity condition as follows. If (JM −1 J T + Jv + b) i > 0 holds for some i, then the rigid bodies are moving away from each other in the direction of the ith constraint, therefore, the ith constraint should be "inactive". However, is the impulse vector, thus i > 0 implies that the ith constraint must be "active". Hence, i > 0 and (JM −1 J T + Jv + b) i > 0 should not hold simultaneously, and this requirement is implemented in the complementarity condition.
By denoting q = Jv + b and A = JM −1 J T and introducing an auxiliary variable w ∈ ℝ m , the LCP (1) can be expressed in a general LCP as follows: It is known that if the constraints are non-degenerate, the Jacobian matrix J is full row rank and the matrix A = JM −1 J T is positive definite (see [5], for example). Throughout this paper, we assume that A is positive definite. Since any positive definite matrix is a P-matrix, A is also a P-matrix and therefore the LCP (3) has a unique solution for an arbitrary vector q.
At the end of this subsection, we should note that the input data b, v, J and M vary in accordance with simulation steps. If we express the time dependence explicitly, they should be b (t) , v (t) , J (t) and M (t) where t is the simulation step. However, in this paper, we mainly focus on solving (1) in each simulation step, therefore, we usually drop the simulation step (t) from b (t) , v (t) , J (t) and M (t) .

Projected Gauss-Seidel method
In the LCP (3) from the rigid-body simulation, the matrix A has a structure such that A = JM −1 J T . The projected Gauss-Seidel (PGS) method [9] is designed to solve more general LCPs (3) in the sense that the assumption for A is only positive definiteness. The PGS method is an iterative method and generates a sequence k ∞ k=0 ⊂ ℝ m . A key property in the PGS method is to decompose A into A = D − L − U such that D is a diagonal matrix, L a strictly lower triangular matrix, and U a strictly upper triangular matrix. Since we assume A is a positive definite matrix, D is invertible. Due to this decomposition, A + q = 0 is equivalent to = D −1 (L + U − q).
Taking this formula and the complementarity condition into consideration, the PGS method computes the next iteration k+1 by the following update formula: Throughout this paper, we use max {a, b} ( min {a, b} ) to denote the element-wise maximum (minimum, respectively) of two vectors a and b . The PGS method continues the update by (4) until the sequence k ∞ k=0 converges enough, or the number of the iterations reaches a certain limit.
In the rigid-body simulation, the initial vector 0 is usually set as a zero vector 0 or the impulse vector obtained in the previous simulation step. Since the initial vector often affects the performance of iterative approaches, the use of the solution from the previous simulation step makes the convergence faster [9]. Such a technique is called warm start.

Accelerated modulus-based Gauss-Seidel method
For solving the general LCP (3), Bai [2] devised the following implicit fixed-point equation that is essential for the modulus-based matrix splitting iteration (MMSI) methods: Here, the matrices M 0 ∈ ℝ m×m and N 0 ∈ ℝ m×m are a splitting pair of A such that A = M 0 − N 0 . ∈ ℝ m×m and ∈ ℝ m×m are two diagonal matrices whose diagonal entries are positive. 1 ∈ ℝ m×m and 2 ∈ ℝ m×m are nonnegative diagonal matrices such that = 1 + 2 . We should emphasize that the variable in (5) is x , and we use |x| to denote the element-wise absolute values of x . The relation between x and the pair of and w in (3) will be discussed in Theorem 1.
By setting 1 = , 2 = O , and = 1 E m with a parameter > 0 in (5), we can obtain a simplified implicit fixed-point equation: Based on (6), the iteration of the MMSI method can be derived as follows: We decompose A into A = D − L − U in the same way as the PGS method. Set = 2 , and let > 0 and > 0 be two parameters. Then, we can derive the update formula of four methods from (7); the modulus-based Jacobi (MJ) by setting M 0 = D in (7) Zheng and Yin [20] established the following theorem to show an equivalence between (8) and LCP(q, A) in (3). Since a detailed proof is not given in [20], we give the proof here.
Proof We first prove (i). Since ( , w) is a solution of LCP(q, A) , ( , w) satisfies the four constraints, A + q = w , w T = 0 , ≥ 0 and w ≥ 0 . The first constraint A + q = w is equivalent to From the rest three constraints and the fact that and are diagonal matrices whose diagonal entries are positive, if x = 1 2 ( −1 − −1 w) , it holds that |x| = 1 2 ( −1 + −1 w) . Therefore, x satisfies and this is equivalent to (8). ) . B y t h e r e l a t i o n s = (|x| + x) and w = (|x| − x) , we obtain A + q = w . Since and are positive diagonal matrices, it is easy to check that and w are nonnegative vectors. Finally, it is also easy to show the element-wise complementarity between and w . ◻ We may use Theorem 1 to establish some iterative methods for solving LCP(q, A) , but we need to set appropriate matrices for the implicit fixed-point equation (8) in actual computations. In particular, the splitting pair of is not unique. By fixing 1 = , 2 = O and = 1 E m , we derive a simplified update equation of (8) as follows: As mentioned in Sect. 2.1, we use E r ∈ ℝ r×r to denote the identity matrix of order r. Based on this equation, Zheng and Yin [20] provided an update formula of the AMMSI methods: When the sequence x k ∞ k=0 converges enough, the AMMSI methods output the impulse vector by using the relation = (|x| + x) = |x|+x . As Mezzadri [11] points out, every AMMSI method can be written in a projection form, and the value of does not play an important role in convergence rates as it just changes the first iteration of the method in its projection form.
By changing the splitting pairs of A , the update formula (11) above yields variant methods; MMSIM ( M 2 = A and N 2 = O ), the accelerated modulus-based Jacobi (AMJ) iteration method ( M 1 = D , N 1 = L + U , M 2 = D − U and N 2 = L ), the accelerated modulus-based SOR (AMSOR) iteration method ( , and the accelerated modulusbased accelerated overrelaxation (AMAOR) iteration method In particular, the update formula of the accelerated modulus-based Gauss-Seidel (AMGS) method in [20] is Then, (12) is equivalent to thus k is a multiple of the positive part of x k . This motivates us to split x k into the positive and negative parts s u c h t h a t (14) is equivalent to Therefore, for computing x k i , we only need the ith component of x k − , which will be denoted as (x k − ) i , since D and are diagonal matrices. This simplifies the computation of x k . Recalling the decomposition of A = D − L − U , we compute x k i for each i = 1, … , m by We can summarize a framework of the AMGS method as follows.

Method 1 (the AMGS method for (3))
Choose a nonnegative vector 0 ∈ ℝ m as an initial vector. Generate the iteration sequence { k } ∞ k=0 by the following procedure:

Convergence theorem
In this section, we focus on the convergence of the accelerated modulus-based Gauss-Seidel (AMGS) method.
Although Zheng and and Yin [20] gave a generic form of a sufficient condition of the convergence, a more detailed discussion of the case that A is a positive definite matrix is somewhat limited: is assumed to be a multiple of the identity matrix ( =̄E m for some ̄> 0 ), so we cannot take = D that is actually used in the numerical experiments in [20].
Since the matrix A is always positive definite in the rigid-body simulation, it is worthwhile to give a sufficient condition of the convergence when A is positive definite while allowing a more generic form of . For example, if we can take = D , we may be able to improve the convergence, as we will show in Sect. 5. Here, ∈ ℝ is a positive constant, and D is the diagonal matrix whose diagonal elements are those of A . We also give an example in Example 1 that is covered by the theorem we will show, but is not covered by theorems in [20].
We use ‖x‖ = √ x T x to denote the Euclidean norm of x ∈ ℝ m . We also use ‖L‖ to denote an arbitrary norm of a matrix L ∈ ℝ m×m that is consistent with the Euclidean vector norm, so that we can employ ‖Lx‖ ≤ ‖L‖‖x‖ . Especially, ‖L‖ 2 represents the spectral norm of a matrix L.
The following theorem covers the case that = D , which will be actually used in numerical experiments of Sect. 5.
For the subsequent discussion, we assume that A is a symmetric positive definite matrix and hence U = L T . Let ̄ = . The matrix ̄ is also a positive diagonal matrix.

Theorem 2 Let
A be a symmetric positive definite matrix, ̄ be a positive diagonal matrix, and A = D − L − L T be a splitting of the matrix A such that D is a positive diagonal matrix and L is a strictly lower triangular matrix. Also, let ‖⋅‖ be an arbitrary submultiplicative matrix norm consistent with the Euclidean vector norm. Then, if the inequality holds, the iteration sequence { k } ∞ k=0 ⊂ ℝ m generated by Method 1 with an arbitrary nonnegative initial vector 0 ∈ ℝ m converges to the unique solution * ∈ ℝ m of LCP(q, A).

Proof
In the AMMSI methods [20] from which the AMGS was derived, we chose 1 = , 2 = O and = 1 E m for the simplified fixed-point equation (10). Let ( * , w * ) ∈ ℝ m × ℝ m be a solution of LCP(q, A) , and let x * = 1 2 ( * − −1 w * ) . From Theorem 1, x * satisfies (10), and the convergence of k ∞ k=0 to * can be guaranteed by that of x k ∞ k=0 to x * , therefore, we discuss the convergence of the sequence x k ∞ k=0 . Subtracting (10) with x * from the update formula (11) holds, we obtain the convergence of {x k } ∞ k=0 by the fixedpoint theorem and we know that x * is the unique point to which the sequence {x k } ∞ k=0 converges. Note that (17) Based on Theorem 2, we obtain the following corollary. Proof By applying the spectral norm to (17), we get From the fact that holds for any invertible matrix P and its singular values i , we have Thus, follows from the assumption and this completes the proof. ◻

Example 1
The following example gives a case that is applicable to Corollary 1, but is not applicable to Theorem 4.1 in [20]: and 2‖L‖ 2 < min − max holds. However, it holds that so (̄ ) + 2 (̄ ) + 2 (̄ ) = 1.2653 … > 1 and therefore the example does not fulfill the assumption of Theorem 1.4 in [20].

Accelerated modulus-based matrix splitting iteration methods for interactive rigid-body simulation
Since the AMGS method is a general method for LCPs, it is possible to simply apply the AMGS method to (3). However, explicit evaluation of A = JM −1 J T is inefficient even though the matrices J and M −1 are sparse. Thus, such a simple application of the AMGS method is not practical. To overcome this inefficiency, we modify the AMGS method so that it does not require the explicit evaluation of the matrix A . In a similar way to the AMGS method, we can also modify the AMSOR method for solving LCPs in the rigid-body simulations. As already pointed out at above, the direct computation of A is unfavorable for real-time simulations. Thus, we should avoid the computations of (15), which involve all the off-diagonal elements of A.
To improve the computational efficiency, we introduce an intermediate variables k+1,i ∈ ℝ m and v k+1,i ∈ ℝ n . In particular, v k+1,i stores information of applied impulse [9,17]. For i = 1, … , m , let where M = M −1 J T . By the definitions of A and q , it holds and this leads to Due to the relation k = 2 x k + , we can compute k+1,i+1 by updating only the ith position of k+1,i ,

Thus, from (19), we obtain
where M * i is the ith column of M .
The following method summarizes the proposed AMGS method for rigid-body simulations.

Method 2 (the proposed AMGS method for rigid-body simulations)
Choose a nonnegative vector 0 ∈ ℝ m as an initial vector. Generate the iteration sequence { k } ∞ k=0 by the following procedure:

satisfies a certain convergence threshold
In Sect. 5, we compare the computational costs of Method 1 and Method 2, which is an optimized version of Method 1 for interactive rigid-body simulations. As the experiments will show, Method 2 achieves considerably lower computational costs than that of Method 1 in all cases, and is suitable for interactive rigid-body simulations.

Linear complementarity problems with lower and upper bounds
In this section, we discuss a method for solving LCPs with lower and upper bounds on , which often occur in the formulations of contact constraints with friction [14,17]. We call such LCPs with lower and upper bounds "Boxed LCPs (BLCPs)". Consider the following BLCP: Here, l ∈ ℝ m and u ∈ ℝ m are the lower and the upper bounds, respectively. Without loss of generality, we assume 0 ≤ l i < u i for each i = 1, … , m.
We define a projection function of to the interval l and u by Since l is a nonnegative vector, so is p BLCP ( ) . With this projection, we can give an AMGS method for Boxed LCPs as follows.

Method 3 (the AMGS method for Boxed LCPS in rigidbody simulations)
Choose a nonnegative vector 0 ∈ ℝ m as an initial vector. Generate the iteration sequence { k } ∞ k=0 by the following procedure:

satisfies a certain convergence threshold
Note that BLCP(q, l, u, A) with l = 0 and u being a very large vector represents the same problem as LCP(q, A) , thus Method 3 can be considered as a generalization of Method 1.

Numerical experiments
In this section, we show numerical results of several 2-dimensional examples to verify the numerical performance of the proposed method. All tests were performed on an Windows 8 computer with Intel Core i7-5500U (2.4 GHz CPU) and 8 GB memory space.
In the numerical experiments, we utilized the warmstart strategy [9], that is, the final solution obtained in a simulation step will be used as the first guess of the solution in the next simulation step. More precisely, let (t),k ∞ k=0 denote the sequence generated for solving LCP(q (t) , A (t) ) at a simulation step t. We set the number of iterations in a single simulation step to 10, that is, (t),10 is used as the first guess (t+1),0 . The initial point of the entire simulation is set as the zero vector ( (0),0 = 0) . To evaluate the accuracy of (t),k , we use a residual function in [2]: For the AMGS methods, we set = 2 , and ̄ is chosen as Pool 1 Figure 1 displays an example with 221 circles in an area of 6 meters wide. All the circles have the same mass (2 kilograms) and the same radius (21 centimeters). The coefficients of friction are set to 0.1, and the coefficients of restitution are set to 0.2. When the coefficient of restitution is 1, collisions are perfectly elastic (i.e., no energy loss), and when the coefficient of restitution is 0, collisions are perfectly inelastic. The gravitational acceleration is set to 9.80665 m/s 2 in the downward direction in the figure.
In Pool 1, the size n in the matrices J ∈ ℝ m×n and M ∈ ℝ n×n is 672, while the size m depend on the simulation step t, since the number of contact points between

Pool 2
The most part of this example is same as Pool 1, but the masses of the circles increases linearly from 1.0 kilograms

Convergence in each simulation step
In each simulation step t, we execute only 10 iterations and we move to the next simulation step t + 1 with the first guess (t+1),0 = (t), 10 . In this subsection, we execute more iterations to compare the convergence of the PGS method and the proposed AMGS methods in each simula-   . 8 The residuals in the entire simulation of Stacking RES( (100),k ) of Pool 2 and Stacking, respectively. We chose the 100th simulation step, since the early steps contained a lot of noise and the warm-start strategy did not work effectively there. From Figs. 3, 4, and 5 , we observe that Method 2 attains better convergence than the PGS method in most cases. Table 1 reports the iteration number k and the computation time in seconds of each method to reach RES( (100),k ) < 10 −4 . Since the standard AMGS method cannot directly handle contacts with frictions, coefficients of friction are set to 0 only for this experiment. Also, the computation time for the 30 circles in Stacking (Fig. 2) was too short to measure (shorter than 0.001 s), therefore we used 150 circles instead of 30 circles. Method 1 computes the coefficient matrix A explicitly, and we observe that Method 1 is the slowest in Table 1

Convergence in entire simulation
In this subsection, we report the computational errors RES( (t),10 ) along with the progress of simulation steps t ≥ 60 . We removed the first 60 steps from the figure, since the early steps contained much noise. In the previous subsection, we observed that Method 2[ 1 D ] is superior to Method 1 and Method 2[̄E m ] in each simulation step. Thus, we use only Method 2[ D ] in this subsection, changing the value of .
In Fig. 6, the horizontal axis is the simulation step t, and the vertical axis is the residual RES( (t),10 ) . From Fig. 6, we observe that Method 2 converges faster than the PGS method for t ≥ 300 . Among different values of , = 0.2 shows the fastest convergence in the figure.
The result of Pool 2 illustrated in Fig. 7 indicates that there are no clear differences of the convergence speed between the PGS method and the AMGS method with various values of , but when = 0.1 , the simulation is unstable during about the first 100 simulation steps.
From Fig. 8 for Stacking, in a similar way to Pool 1, we can again observe that the AMGS method with the smaller gives the faster convergence, and the AMGS method with = 0.6 still converges faster than the PGS method.
Finally, Table 2 shows the entire computation time for 1,000 simulation steps, with 200 iterations for each step, that is, we computed (1000),200 . The entire computation time is shorter in the proposed method than in the PGS method. As mentioned in Table 1, the convergence rate is better in the proposed method (Method 2[ D ]) than in the PGS method, in other words, the proposed method simulates the circles more accurately than the PGS method. Therefore, the proposed method has the advantages for interactive rigid-body simulations.

Conclusion
We presented a numerical method based on the AMGS method for interactive rigid-body simulations. We established the convergence theorem of the AMGS method for the case the matrix A is positive definite and = D with > 0 . This case was examined in the numerical experiments, and we observed that the proposed method attained the better accuracy than the PGS method and the computation time of the proposed method was shorter than that of a simple application of the AMGS method.
In practical cases, however, determining a proper value of is not simple. The numerical results showed that a smaller value of gave a better convergence, but it is also shown that a too small value of results in an unstable behavior. An approach that adaptively determines the value of may resolve this problem, and we leave a discussion on such an approach as a future task of this paper. Further numerical experiments in 3-dimensional spaces that take frictions into consideration will be another topic of our future studies. Related to computing frictions correctly, applying the proposed method to more general form of LCPs, such as NCPs and HLCPs, will also be an interesting extension of this paper.
Funding This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

Compliance with Ethical Standards
Conflict of interest On behalf of all authors, the corresponding author states that there is no conflict of interest.
Availability of data and material All data generated or analysed during this study are included in the article.

Code Availability Not applicable
Disclosure of potential conflicts of interest On behalf of all authors, the corresponding author states that there is no conflict of interest.

Informed consent Not applicable
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.