Symmetric and Asymmetric Multiple Impulsive Constraints Without Friction and Their Characterization

We present two meaningful and effective non-ideal constitutive characterizations for a multiple impulsive constraints S\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {S}}$$\end{document} comprising a finite number of non-ideal frictionless constraints of codimension 1, described in the geometric setup given by the space–time bundle M\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {M}}$$\end{document} of a mechanical system in contact/impact with S\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {S}}$$\end{document}. Thanks to the geometric structures associated to the elements of S\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {S}}$$\end{document}, we introduce a symmetric characterization, that does not distinguish the elements forming S\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {S}}$$\end{document} as regards mechanical behavior, and an asymmetric one that makes this distinction. Both the characterizations provide a generalization of the characterization of ideal multiple constraints presented in Pasquero (Q Appl Math 76(3):547–576, 2018). The iterative nature of these characterizations allows the introduction of two algorithms determining the right velocity of the system in case of single or multiple contact/impact with symmetric or asymmetric constraints S\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {S}}$$\end{document}, once the elements forming S\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {S}}$$\end{document} and the left velocity of the system are known. We show the effectiveness of the two possible choices with explicit implementations of these algorithms in two significant examples: a simplified Newton’s cradle system for the symmetric characterization and a disk in multiple contact/impact with two walls of a corner for the asymmetric one.

and more specifically in classical impulsive mechanics since the early 2000 (see e.g., Pasquero 2005).
One of the best assets of jet bundle theory applied to classical mechanics consists in its natural framing of time-dependent description and frame independent description of the Lagrangian representation of mechanical phenomena. Frame independence is particularly important in impulsive mechanics, where the concept of velocity, intrinsically dependent on the choice of a frame of reference, plays a crucial role in the laws that govern the motion both for free and constrained systems.
The jet bundle geometric context are especially suitable to frame the so-called event-driven algebraic approach to impulsive systems, where every impact happening in the time evolution of the system is considered as an isolated event with the system in a fixed configuration in contact with the constraints and where the kinematic state of the system after the impact can be expressed as function of the kinematic state before the impact. In particular, the jet bundle context allows a very simple definition of constitutive characterization of impulsive and/or unilateral constraints, that is, roughly speaking, the rule to which the reactive impulse is subject in order to restore determinism. On the understanding that the only compulsory condition for a characterization is to restore determinism of classical mechanics, the choice of a characterization can be based on several different requirements: for instance geometric properties of the constraint itself, frame independence, ideality, extremality, generalization of known cases, physical meaningfulness, agreement with experimental data and so on. The following different cases of impulsive constraints and characterizations have found a satisfactory description in the geometric context of jet bundle theory: • Ideal single positional constraints of codimension 1 and of codimension greater than 1 (see Pasquero 2005); • Ideal constraints of mixed nature (positional with additional permanent and/or instantaneous kinetic constraints; see Pasquero 2006); • Non-ideal single constraints of mixed nature (this case can be divided into nonideal constraints without or with friction; see Pasquero 2008); • Ideal multiple positional constraints formed by 2 ≤ r < n constraints of codimension 1 (see ). The relationships between the cases listed above are illustrated by the following diagram that shows, starting from the "progenitor" case of ideal single positional constraint of codimension 1, the increasing refinement of the modeling and the essential distinction between single and multiple constraints. Multiple contact/impact, that is the situation when a mechanical system is simultaneously in contact with two o more constraints and impacts with one or more of them, occurs for many mechanical systems, from scholarly cases such as billiard balls to more up to date and applicable cases such as granular matter and kinematic chains.
The aim of this paper is to present one step further in the left column of the diagram above by defining and analyzing two possible constitutive characterizations of nonideal frictionless multiple positional constraints in symmetric and asymmetric cases. In the symmetric case, the characterization consists in a direct generalization of the ideal one by introducing a suitable global coefficient of restitution (COR) applicable to each element of the multiple constraint. In the asymmetric case, the characterization differentiates the single constraints by considering a (possibly) different COR for each element of the multiple constraint. Of course, the asymmetric characterization comes down to the symmetric one when the CORs of the single constraints are all the same.
The geometric context to which we refer in this paper is the same of Pasquero (2018), fit to manage multiple positional constraints formed by 2 ≤ r < n constraints of codimension 1. In order to avoid the technicalities of multiple indexes, and without significant loss of generality, we will consider r = 2, then modeling simultaneous contact/impact of the mechanical system with two positional constraints of codimension 1. Moreover, we restrict our attention to systems subject only to positional (and not kinetic) constraints.
With the intention of getting the paper closer to self-consistency, and also in order to fix notation, a very concise description of the geometry of the Lagrangian setup of constrained impulsive mechanical systems subject to multiple positional constraints is presented in Sect. 2. It is a short version of that presented in . The reader interested in a wider survey on the argument can also refer to .
The main known ideas about the concept of impulsive constitutive characterization and the application to single constraints and ideal multiple constraints are presented in Sects. 3.1 and 3.2, respectively. The characterizations of symmetric and asymmetric non-ideal frictionless multiple constraints, that constitute the innovative part of the paper, are presented in Sects. 3.3 and 3.4. The general procedure applicable when r > 2 is sketched in Sect. 3.5.
The iterative nature of the procedure giving the characterizations allows the implementation of algorithms effectively applicable to study impacting mechanical systems: Sect. 4 presents two examples of these applications. In Sect. 4.1, we illustrate the behavior of the symmetric characterization applied to a simplified version of the Newton's cradle, formed by three equal balls. This is a case of "multiple contact with repeated single impacts in a symmetric constraint" (different to the case "multiple impact in a symmetric constraint," an example of which is analyzed in Fassino and Pasquero (2019)). In Sect. 4.2, we illustrate the behavior of the asymmetric characterization applied to a "multiple impact in an asymmetric constraint," analyzing the multiple impact of a disk in a corner formed by two walls having different CORs.
The list of references is based on the minimality criterion of making the paper self-consistent and then is focussed on the specific approach given by the jet bundle framework. Different choices, even if restricted to the works pertaining only multiple impacts, could draw away the attention of the reader from the peculiar approach of this paper and moreover could give rise to reasonable criticisms and concern (Saracco 2020). However, the reader interested in classical approach to contact/impact mechanics can find foundations, main results and a huge bibliography in Johnson (1985), Stronge (2000), Brogliato (1996). A different approach to frictionless multiple impacts can be found in Liu et al. (2008Liu et al. ( , 2009. Different but still heavily based on geometric techniques approaches to unilateral and impulsive constraints can be found in Ibort and De Leon (1998), Cortés and Vinogradov (2006).

Geometry of the Lagrangian Setup
In this part, we concisely describe the Lagrangian geometric setup suited to study impulsive mechanics of systems subject to positional constraints, and we recall the main definitions and results about the characterization of impulsive constraints in this context. With few exceptions, the material is not new (for a more detailed discussion we refer to  and the references within).

Free Systems
The configuration space-time of a mechanical system with a finite number n of degrees of freedom is a fiber bundle π t : M → E where M is a (n + 1)-dimensional differentiable manifold, E is the 1-dimensional euclidean space representing the time axis, and π t is the projection implementing the absolute time axiom. We usually refer M to fibred coordinates (t, x 1 , . . . , x n ), where t is a global euclidean coordinate on E. A motion of the system is a section γ : E → M, locally represented by a map γ = γ (t) such that t (t, x 1 (t), . . . , x n (t)). The (2n + 1)-dimensional vector bundle π : V (M) → M of the vertical (with respect to π t ) vectors that are tangent to the fiber of M is the subbundle of the tangent bundle π : T (M) → M given by the space-like vectors. The elements of V (M) have local representation V = V 1 ∂ ∂ x 1 + · · · + V n ∂ ∂ x n . The (2n + 1)-dimensional affine bundle π : J 1 (M) → M of the vectors that are tangent to any possible motion of the system in any point is the subbundle of the tangent bundle π : T (M) → M given by the time-like vectors, also called absolute velocities.
where the (mass) matrix g i j takes into account of the massive properties of the system.
is the kinetic energy of the system with respect to the frame h.
An impulse acting on the system is an element I ∈ V (M). The mechanical law governing the impulsive phenomenon is simply the action + : given a left velocity p L of the system and an impulse I such that π(p L ) = π(I) ∈ M, we can determine the right velocity shows the independence by the frame of reference of the jump of velocities given by I.

Positional Constraints
A positional constraint S is a (s+1)-dimensional fibred subbundle i : S → M without boundary. The bundle S determines the following additional geometric objects and structures relative to the system: • The contact condition of the system with the constraint. The system is in contact with the constraint if its configuration x ∈ i(S) ⊂ M. Since we adopt the eventdriven approach; later on the system will be always considered in contact with the constraint; representing the bundles of time-like and space-like vectors of the system when the system is in contact with the constraint S; Some comments are in order to apply these structures in the context of impulsive mechanics.
Heuristically, the pull-back bundle (i * ) * (J 1 (M)) is formed by the absolute velocities of the system p ∈ J 1 (M) applied in points of S but not necessarily tangent to S: then (i * ) * (J 1 (M)) is the natural geometric framework fit to analyze the behavior of S viewed as a unilateral positional constraint. Since the system is always supposed in contact with the constraint, focussing mainly on the fibers of the bundles, later on we make the identifications i * (J 1 (S)) J 1 (S), i * (V (S)) V (S), and the identi- V (M) (even though these last identify bundles that are structurally different since they have different base manifolds).
Given an absolute velocity p ∈ J 1 (M) of the system in contact with S, the vertical vector V ⊥ S (p) = P ⊥ J 1 (S) (p) assumes the "frame invariant" meaning of orthogonal component of the absolute velocity p with respect to S. We recall, however, that the concept of tangent component of p with respect to S has not a frame invariant meaning (see Pasquero 2005).
A positional constraint S is called unilateral if two sets λ(S) ⊂ J 1 (M) and ρ(S) ⊂ J 1 (M) are assigned so that J 1 (M) can be written as the disjoint union (1) If codim(S) = 1, then S has a natural structure of unilateral constraint. Since Then, for instance, we can define

Multiple Positional Constraints
A multiple positional constraint S is the assignment of two (or possibly more) positional constraints S 1 , S 2 satisfying The system has a (multiple) contact with S if its space-time configuration x ∈ S. Due to condition (iii) above, a multiple positional constraint can also be considered as unilateral. For instance, labeling for simplicity with the indexes 1, 2, 12 the objects relative to S 1 , S 2 , S, respectively, we can set: For every p ∈ J 1 (M), we have then three different orthogonal velocities: The following result clarifies the relation between these velocities.

Proposition 1 For every
Then, due to the independence condition (iii) and the dimensions of the spaces, we Since P 12 (p) = P 12 (P 1 (p)) = P 12 (P 2 (p)), we have that Therefore, by orthogonality, we have λ = μ = 1 and so the first statement. About the second statement, we have that: The absolute velocity p determines a multiple impact if p ∈ λ(S 1 ) and p ∈ λ(S 2 ) (and of course it determines a single impact with multiple contact if the space-time configuration π(p) ∈ S and p ∈ λ(S 1 ) "exor" p ∈ λ(S 2 )).

Constitutive Characterizations of Impulsive Constraints
In the frame independent context of Lagrangian systems, a characterization for an impulsive constraint (not necessarily positional and/or of codimension 1) is a map I : J 1 (M) → V(M) assigning to every absolute (left) velocity p L ∈ J 1 (M) of the system in contact with the constraint a uniquely determined reactive impulse, that is a space-like vector I(p L ) ∈ V (M). Once I(p L ) is known, an absolute (right) velocity p R ∈ J 1 (M) is determined by the rule p R = p L + I(p L ) (see e.g., . Later on the constraints will be always considered as unilateral.

Single Constraints
The simplest case of positional constraint S with codim(S) = 1 is emblematic: The only requirements of frame independence and ideality, in the form of conservation of kinetic energy for all the rest frames of S, univocally determine the characterization as the "reflection rule" I(p L ) = −2 V ⊥ S (p L ) (see Pasquero 2005). An analogous characterization turns out to be significant also in more complex cases, such as ideal constraints of codim(S) > 1 or in presence of ideal kinetic constraints of both permanent or instantaneous kind (see Pasquero 2005Pasquero , 2006. The introduction of a (Newtonian) COR ε ∈ [0, 1] allows to generalize the "reflection rule" also in case of non-ideal frictionless constraints, in the form Pasquero 2008). Of course, the condition ε = 1 returns the case of ideal constraint, while the condition ε = 0 gives in this case the totally inelastic impact, with the annihilation of the normal (with respect to S) component of the impact velocity and the maximum possible loss of kinetic energy of the system (with respect to every frame h S ∈ H S ). Note, however, that this does not imply the rest of the system after the impact, since this concept has a clear meaning only if a frame of reference is fixed.
To give a physical meaning to the unilaterality of the constraint S, the characterizations above must be detailed in the form (with ε ∈ [0, 1]) This assignment naturally gives an outgoing right velocity p R starting from an incoming left velocity p L .

Multiple Ideal Constraints
Due to the presence of three different orthogonal velocities, the same ideas are not analogously appropriate to model a multiple constraint S = S 1 ∩ S 2 . Once again labeling for simplicity with the indexes 1, 2, 12 the objects relative to S 1 , S 2 , S, respectively, the simplest rule I(p L ) = −2 V ⊥ 12 (p L ) does not take into account the multiple nature of the constraint, while a rule combining V ⊥ 1 (p L ) and V ⊥ 2 (p L ) must take into account several details, for instance: • The "incoming" nature of p L , that could be incoming for S 1 alone, or for S 2 alone (in case of single impact), or for both (in case of multiple impact); • The symmetric or asymmetric behavior of the constraints S 1 and S 2 (that is if the two constraints have or not the same COR); • The verification of the outgoing nature of p R .
An ideal characterization for multiple constraints p ∈ λ(S 1 ) and p ∈ λ(S 2 ). See Pasquero 2018) can be assigned on the basis of the rule where the coefficient β is Note that, since by Prop.
(1) we have β = 1 if V ⊥ 1 (p L ) = 0 "exor" V ⊥ 2 (p L ) = 0, the rule restores the usual reflection law in case of single constraint. Moreover, the rule restores also the usual characterization I(p L ) = −2 V ⊥ 12 (p L ) of single constraints of codimension greater than 1 (see Pasquero 2005) in case of orthogonality V ⊥ 1 (p L ), V ⊥ 2 (p L ) = 0. Thirdly, note that the coefficient β has a "kinematic" dependence on p L , and not only a "geometric" dependence on U ⊥ 1 , U ⊥ 2 . However, even in case of single impact with multiple contact, the absolute velocity p L + I(p L ) could be not an exit velocity, so that, in order to have an outgoing p R , the evaluation of the reactive impulse could (possibly) be iterated starting from the new entrance velocity p L + I(p L ) (see .

Symmetric Multiple Non-ideal Frictionless Constraints
The ideal characterization (5) allows, once a COR ε is assigned, a natural generalization in the form Such a rule models the symmetric non-ideal behavior of the constraint S: In fact the rule, handling symmetrically the two constraints, restores the usual non-ideal frictionless characterization in case of single impact for both the single impacts with S 1 and S 2 .
Once again the absolute velocity p L + I(p L ) could be not an exit velocity, so that the evaluation of the reactive impulse could (possibly) be iterated.
In the next section, we will analyze the well-known system given by a simplified Newton's cradle formed by three equal balls. This is an example of system subject to multiple contact and (iterated) single symmetric impact. The application of this rule in the case of a disk impacting with two equal walls forming a corner, that is a case of a system subject to multiple contact and multiple symmetric impact, is analyzed in details in Fassino and Pasquero (2019).

Asymmetric Multiple Non-ideal Frictionless Constraints
Both the characterizations (5, 6) can be generalized, once two possibly different CORs ε 1 , ε 2 for S 1 , S 2 , respectively, are assigned in the form Such a rule models the asymmetric non-ideal behavior of the constraint S. The rule restores the usual non-ideal frictionless characterization in case of single impact with CORs ε i for the single impacts with S i . Moreover, it restores the previous characterization (6) in the case ε 1 = ε 2 and the ideal one (5) for ε 1 = ε 2 = 1. Once more, however, the absolute velocity p L + I(p L ) thus obtained could be not an exit velocity and the evaluation of the reactive impulse could (possibly) be iterated. The application of this rule to the case of a disk impacting with two walls (having different CORs) forming a corner is analyzed in details in the next section.

Remarks on the Characterizations
The (iterative) line of reasoning traced above to assign the right velocity p R starting from a left velocity p L can be formalized also for a multiple constraint S = r i=1 S i formed by 2 ≤ r < n single constraints with the following step by step procedure. Given a left velocity p 0 L ∈ J 1 (M) and the CORs ε 1 , . . . , ε r ∈ [0, 1]: a) Select the set of "contact" indexes {i 1 , . . . , i k } ⊂ {1, . . . , r } such that the configuration π(p 0 L ) ∈ M is such that π(p 0 L ) ∈ S i m ∀ m = 1, . . . , k. If {i 1 , . . . , i k } = ∅ then let p R := p 0 L (the system is not in contact with the constraint). Otherwise b) Select the set of "impact" indexes { j 1 , . . . , j h } ⊂ {i 1 , . . . , i k } such that the orthogonal velocities V ⊥ j q (p 0 L ) are such that (p 0 L ∈ λ(S j q ) ∀ q = 1, . . . , h); c) If h ≥ 2, calculate the "global" orthogonal velocity V ⊥ j 1 ... j h and the corresponding coefficient

Example 1. The Simplified Newton Cradle
Three equal disks of radius R and unitary mass can move in a plane. Labeling the disks with the numbers 1, 2, 3, the space-time configuration is described by 10 coordinates (t, x i , y i , ϑ i ), i = 1, . . . , 3 where (x i , y i ) are the coordinates of the center of the ith disk and ϑ i is its orientation. The set of constraints is given (with obvious notation) by the functions We consider the particular case when disks 2 and 3 are at rest and in contact, while disk 1 moves in the straight line direction given by the centers of 2 and 3 and collides with disk 2 (see the upper part of Fig. 1), so that at the impact instant we have x 2 − x 1 = x 3 − x 2 = 2R and y j − y i = 0, ∀ i, j = 1, . . . , 3. Then the positional constraints involved in the impact are S (1,2) and S (2,3) , and the system is subject to multiple contact.
Starting with left velocity p 0 it is easy to show that, due to the nature of the constraints involved in the impact, all the velocities evaluated by the theoretical algorithm do not have components in the ∂ ∂ y i directions. Moreover, due to the absence of friction between the disks, a possible initial spin ω 0 ∂ ∂ϑ 1 of disk 1 (or even disks 2 and 3) does not have significant effect on the impact. (It will be found unchanged in the final result of the algorithm.) Then the system has only (t, x 1 , x 2 , x 3 ) as significant coordinates, the termsẋ i of the velocity (p = ∂ ∂t +ẋ i ∂ ∂ x i ) represents the horizontal component of the velocity of the ith disk, and the impact conditions could be easily checked in the formẋ r +1 −ẋ r < 0 for some r = 1, 2. Anyway, recalling (2), we have so that the velocity p 0 L is an incoming velocity only for the constraints S (1,2) , and the system is subject to a single impact. The standard calculation of V ⊥ (1,2) (p L ) gives Applying the (first step of the iterative) method described in the previous section, the "new" left velocity is and we have to determine the nature of the new left velocity p 1 L with respect to the set of constraints. We have: so that the new left velocity is an outgoing velocity for S (1,2) and an incoming velocity for S (2,3) . Once again, the system is subject to a single impact. The standard calculation from which we obtain the "new" left velocity and we have to restart. We have: Then, p 2 L in an exit velocity for the system if and only if ε = 1. Otherwise, p 2 L is an outgoing velocity for S (2,3) , but it is once again an incoming velocity for S (1,2) . This "second impact" between disks 1 and 2 is a significant difference with respect to the behavior of the ideal Newton's cradle (see .
Further iterations of the algorithm then present different possibilities depending on the value of ε: the next step gives that is an exit velocity for the system if and only if ε ≥ 3−2 √ 2 0, 17157; otherwise, we have a "second impact" between disks 2 and 3 and the procedure must be further iterated.
With the intent to exit from a completely theoretical presentation (similar to that of Pasquero (2018)) of the characterization (6) applied to the Newton's cradle with three equal balls with initial velocity p 0 L = ∂ ∂t + v 0 ∂ ∂ x 1 , the algorithm for the determination of the exit velocity of the system has been implemented using the software PariGP (The PARI Group (2019). The PariGP implementation of the algorithm for the simplified Newton's cradle is available as Electronic Supplementary Material). In order to avoid to delve into the technicalities of the robustness of the algorithm (such as in Fassino and Pasquero 2019), the termination was ensured by fixing a maximum number N = 10 4 of iterations and a minimum threshold σ = 10 −9 v 0 for the differences of the velocities of the three disks. Moreover, a test checking the possible presence of multiple impacts was included, always giving negative result for these initial data.
Since the term v 0 can be factored in the output velocity, Table 1 lists the output of the algorithm (the last rows rounded to the tenth significant digit) for the fixed value v 0 = 100 and for different CORs ε. The data show a meaningful physical behavior of the system. For instance, the number of iterations increases for decreasing CORs, and the exit velocities of the three disks are increasingly similar for decreasing CORs. In particular, the acronym AEEV in the last three rows stands for "Almost Equal Exit Velocities" followed by the number of iterations before checking the almost equalities of the velocities of the disks (note in fact that the impact conditionsẋ r +1 −ẋ r < 0 for some r = 1, 2 is still verified).
The lower part of Fig. 1 illustrates the behavior of the system after the impact for the three cases ε = 0.9, ε = 0.5, ε = 0.1.

Example 2. The Disk in the Corner
A rigid disk of radius R and unitary mass moves in the part of a horizontal plane delimited by two walls S 1 , S 2 forming an angle 2α ∈ (0, π). The system can be described using local coordinates (t, x, y, ϑ) where x, y are the coordinates of the center of the disk and ϑ is the orientation of the disk. If k = tan α > 0, the walls can be described by the Cartesian relations S 1 : kx − y = 0, S 2 : kx + y = 0. We assume that the disk is in contact with both the walls, so that (x, y) = (− R sin α , 0) (see Fig. 2). Given an absolute velocity p = ∂ ∂t +ẋ ∂ ∂ x +ẏ ∂ ∂ y +θ ∂ ∂ϑ , we have with corresponding coefficient β = 1 + k 2 2 k 2ẋ 2 +ẏ 2 k 4ẋ 2 +ẏ 2 . In analogy with Fassino and Pasquero (2019), we divide the space of absolute velocities of the system into four different zones Z 0 , Z 1 , Z 2 , Z 12 with p ∈ Z i if it determines an impact with S i . Taking into account (3), we have Once the CORs ε 1 , ε 2 of S 1 , S 2 , respectively, and an initial velocity Once again with the intent to exit from a completely theoretical presentation, the algorithm for the determination of the exit velocity of the disk from the angle has been implemented using the software PariGP. (The PariGP implementation of the algorithm for the disk in the corner is available as Electronic Supplementary Material). Once again, the termination was ensured by fixing a maximum number N = 10 4 of iterations and a minimum threshold σ = 10 −9 v 0 (with v 0 = 100) for the norm of exit velocity of the disk. Moreover, since the coefficient β changes with the iterations of the algorithm, a test checking the possible excessive downsizing of the denominator of β was included.
The analysis of the output given by the implemented algorithm for different angles α, different impact angles ϕ and for different CORs ε i shows that • The results in case of ε 1 = ε 2 coincide with those presented in Fassino and Pasquero (2019) (up to the different rounding of the implementations of the algorithm); • The results present the correct symmetry with respect to the horizontal axis; • The results show natural behaviors of the system, such as decreasing norm of the final velocity for decreasing values of the CORs or an increasing number of iterations (i.e., "rebounds") as the angle α between the walls becomes narrower.