Synchronization of two coupled pendula in absence of escapement Federico Talamucci arXiv:1602.06382v1 [math-ph] 20 Feb 2016 DIMAI, Dipartimento di Matematica e Informatica “Ulisse Dini”, Università degli Studi di Firenze, Italy tel. +39 055 2751432, fax +39 055 2751452 e-mail: federico.talamucci@math.unifi.it 2010 Mathematics Subject Classification: 34C15, 34L15, 70E55 Keywords: Synchronization, coupled pendula, characteristic equation, eigenvalue localization. Abstract. A model of two oscillating pendula placed on a mobile support is studied. Once an overall scheme of equations, under general assumptions, is formulated via the Lagrangian equations of motion, the specific case of absence of escapement is examined. The mechanical models consists of two coupled pendula both oscillating on a moving board attached to a spring. The final result performs a selection among the peculiar parameters of the physical process (lenghts, ratio of masses, friction and damping coefficients, stiffness of the spring) which provide a tendency to synchronization. 1 Introduction Systems of coupled oscillations are largely studied on account of their wide possibility of application in many significant branches (mechanics, medical and byological sciences, ...). The corresponding mathematical problem is in no way easy to handle when all the effects are overlapped: here, we propose a basic situation which will be discussed from the mathematical point of view. The main question we deal with is the feasibility of in–phase or antiphase synchronization when no external forces (escapement) forcing the free oscillation are contemplated. We mainly take care of the mathematical path drawn by the equations of motion, aiming at developing the analytical scheme, even if in a simplified situation. We first formulate the mathematical model by allowing very general features of the mechanical phenomenon, admitting different sizes of masses, lenghts of the pendula and including escapement conditioning the oscillations. This will supply a ground in order to make a brief comparison with some significant models proposed in literature. Rather than obtaining information directly via a numerical simulation approach, we rather aim for via an analytical method of locating the eigenvalues linked with the damping of the system. The advantage is the prospect of recognizing some ranges of the parameters entering the phenomenon which predispose the sytem to synchronization. On the other hand, some expected results are confirmed, as the unattainability of the in–phase synchronization in absence of escapement. The first step, introduced in the next Section, is the mathematical formulation of the model, which will be achieved by the Lagrangian’s method of deducing the equations of motion. 2 The mathematical model The system we are going to study can be realized either by device I or device II shown in Figure 1: as for I, the apparatus consists of two pendula whose pivots (points A and B) are fasten on a horizontal 1 and homogeneous beam with mass m0 and P0 as centre of mass. The beam is placed on a pair of rollers of radius R. The massive bobs (of mass m1 and m2 ) are suspended at the extremities P1 and P2 of the massless rods (of lenghts ℓ1 < R and ℓ2 < R) and can oscillate on the vertical plane containing the beam. The distances d1 , d2 of P1 and P2 from P0 are larger than ℓ1 and ℓ2 , in order to avoid hits. A spring of stiffness k and whose mass and lenght at rest are negligible, is attached at one of the extremities of the beam. Device I Device II Figure 1: The board supporting the two pendula leans against a couple of rollers (device I) or moves back and forth on a horizontal support (device II). In both cases, a spring is attached to one of the exteremities of the board. In device II the pendula are suspended by means of a T–shaped support, with negligible mass and fixed on the beam which oscillates at a lower height. If the physical quantities m0 , m1 , m2 , ℓ1 , ℓ2 are the same in I and in II, it is immediate to realize that the mathematical problem is identical (actually the lenghts h, R, d1 and d2 do not enter the equations of motion). On the other hand, not even from the dynamical point of view the two systems are different: actually, active forces are the same and, assuming that the rolling friction in apparatus I is proportional to Ψ̇, with Ψ rotation angle of the roll, then it is proportional to ẋ, with x abscissa of P0 (namely x − 2RΨ = constant), just like in apparatus II. 2.1 The equations of motion In this first part the equations of motion are achieved by using a Lagrangian approach. The cartesian coordinate system is fixed in order that the mechanism is contained in the vertical plane y = 0 and the beam swings along the x–axis, the z–axis is upward–vertically directed, the origin O corresponds to the fixed extremity of the spring. We choose as lagrangian coordinates q = (x, θ1 , θ2 ), where θ1 and θ2 are the amplitudes of oscillation with respect to the downward–vertical direction and x is the abscissa of P0 : the representative vector of the discrete system (P0 , P1 , P2 ) in terms of them is, for device I, X(q) = (x, 0, 2R, x − d1 + ℓ1 sin θ1 , 0, 2R − ℓ1 cos θ1 , x + d2 + ℓ2 sin θ2 , 0, 2R − ℓ2 cos θ2 ) 2 and 0 at the third position, h replacing 2R for device II. In both cases, the Lagrangian function L = T +U , where U is the potential of the elastic force and of gravity, is L(q, q̇) = where m = 2 P  2    X 1 1 1 mẋ2 + mj ℓ2j θ̇j2 + mj ℓj ẋθ̇j + g cos θj − kx2 2 2 2 j=1 (1) mj is the total mass. j=0 As for the friction forces, if the damping is formulated as ΦPi = −βi Ṗi , i = 0, 1, 2, the lagrangian components are   2     X βj ℓj θ̇j cos θj , −β1 ℓ1 ẋ cos θ1 + ℓ1 θ̇1 , −β2 ℓ2 ẋ cos θ2 + ℓ2 θ̇2  Φ(q) = −β ẋ − (2) j=1 with β = 2 P βj . If, on the contrary, one assumes that the pendula run into damping only along the j=0 rotational direction eθi = (cos θi , 0, sin θi ), i = 1, 2, then the generalized friction force reduces to   2 X βj ℓj θ̇j cos θj , −β1 ℓ21 θ̇1 , −β2 ℓ22 θ̇2 ,  . Φ(q) = −β0 ẋ − (3) j=1 The mechanism of escapement can be modelled by introducing a moment in the direction j (i. e. orthogonal to the plane of the system), exerting a force fi (Pi , Ṗi , t)eθi on each Pi , i = 1, 2. The corresponding generalized force is F(q) = (f1 ℓ1 cos θ1 + f2 ℓ2 cos θ2 , ℓ1 f1 , ℓ2 f2 ) . (4) d (∇q̇ L) − ∇q L = F(q) + Φ(q) are dt ..    2 2 P P .. mx+ ℓj mj θj cos θj − θ̇j2 sin θj + kx = −β ẋ + ℓj fj − βj θ̇j cos θj , j=1 j=1     .. .. mi ℓi ℓi θi + x cos θi + g sin θi = ℓi fi − βi ẋ cos θi + ℓi θ̇i i = 1, 2 Assuming (2), the equations of motion    2.2   (5) Comparison with some models In reviewing very briefly the mathematical formulation of some models existing in literature, we have the specific intention of remarking that (1) if (3) is accepted to hold, then β0 replaces β in the first equation and the terms −ℓi βi ẋ cos θi , i = 1, 2, 2 P have to be omitted in the second and third equation. However, the term − ℓj βj θ̇j cos θj of the j=1 first equation is present anyway; (2) even if the escapement mechanism operates on the pendula vie the force (4), terms containing f1 and f2 are present in any case in the first equation in (5) concerning x and they affect the motion of the beam. In [2], where the apparatus I is tested, the escapement is formulated by means of a step function, depending on the amplitude of a threshold angle. In [1] the apparatus II is subject to the inversion of direction of the angular velocity (escapement mechanism) at a critical value of the angle. Also in [11] the mathematical problem is formulated for two driven pendula (although the description of the experimental setup refers to a couple of metronomes), but the function (4) is expressed via a continous function. The cited models are undoubtedly significant and useful for the exhibited experimental and numerical results: nevertheless, we remark the lack of the terms in (5), first equation, pointed out in (1), (2) just above. This aspect should not be unimportant, mainly when analytical results are pursued, as in our investigation. 3 An experimental device which differs from I and II described at the beginning of the Section consists in placing two masses at each of the pivotal points and let them oscillate horizontally, by means of a spring connecting the two points (hence the beam is removed): such as interaction mechanism is studied in [3] and the analytical problem is essentially the same as (5), since the centre of mass of the attachement points solves the first of such system, with k = 0. Finally, in the system studied in [5] and [10] the two pendula are coupled by a spring connecting some intermediate points (Q1 and Q2 ) of the two sticks supporting the weights (Kumamoto coupled pendula). The spring–coupled pendula are proposed also as a basic model for the neutrino oscillation. Calling d the constant distance between the pivotal points and presuming that the lenght at rest of 1 1 the spring is d, we write the Lagrangian function as m1 ℓ21 θ̇12 + m2 θ̇22 + m1 gℓ1 cos θ1 + m2 ℓ2 g cos θ2 + 2 2 1 2 k (|Q1 (θ1 ) − Q2 (θ2 )| − d) where 2
1/2 |Q1 − Q2 | = d2 + ℓ20,1 + ℓ20,2 + 2d(ℓ0,1 sin θ1 − ℓ0,2 sin θ2 ) − 2ℓ0,1 ℓ0,2 cos(θ1 − θ2 ) being ℓ0,1 , ℓ0,2 the distances between the pivots and the intermediate points. The problem is here two–dimensional and the corresponding equations of motion for the two angles, by assuming (3) and ℓ1 f1 = u(t), f2 = 0 (external torque acting only on one pendulum, see [10]) are   .. d 2 = u(t) − ℓ21 β1 θ̇12 , m1 ℓ1 θ1 + m1 gℓ1 sin θ1 + ℓ1,0 k[d cos θ1 + ℓ2,0 sin(θ1 − θ2 )] 1 − |Q1 − Q2 |   .. d m2 ℓ22 θ2 + m2 gℓ2 sin θ2 − ℓ2,0 k[d cos θ2 + ℓ1,0 sin(θ1 − θ2 )] 1 − = −ℓ22 β2 θ̇22 |Q1 − Q2 | which slightly differ from the dynamics formulated in [10], where ℓ0,1 = ℓ0,2 and d is neglected. 2.3 The mathematical problem for σ and δ The main purpose is to investigate the existence of solutions of system (5) such that one of the the quantities σ = θ1 + θ2 , δ = θ1 − θ2 (6) tends to zero for t → +∞. If δ → 0 [respectively σ → 0], then the system will proceed to in–phase synchronization [resp. antiphase synchronization].   x Moving toward the more expressive coordinates (6) and setting y =  σ , one has y = Lq, where δ   1 0 1 1 0 1  −1 is the linear change of coordinates. Thus, writing again the equations of motion by   
d d  ∇ẏ Le − ∇y Le = L−1 (∇q̇ L) − ∇q L , F(y) + Φ(y) = L−1 F(q) + Φ(q) with taking account of dt dt ! 1 0 0 e ẏ) = L(L−1 y, L−1 ẏ), L−1 = 0 1/2 1/2 L(y, , we attain L= 0 0 0 1/2 −1/2  .. 1 .. .. + − − + + − σ − +   , Bm )(σ̇ 2 + δ̇ 2 ) + Ψ2 (Bm , Bm )σ̇ δ̇ − Kx+ m x +Ψ1 (Bm , −Bm ) +Ψ1 (Bm , −Bm ) δ = Ψ2 (Bm   2   + − − +    +Ψ1 (f1 + f2 , f2 − f1 ) − β ẋ − Ψ1 (Bβ , −Bβ )σ̇ − Ψ1 (Bβ , −Bβ )δ̇, .. .. .. − + − + − x   σ +A− = −gΨ2 (Bm , Bm ) + f + − Ψ1 (Bβ+ , −Bβ− )ẋ − A+ Ψ1 (Bm , −Bm ) +A+ δ m m  β σ̇ − Aβ δ̇,      .. .. ..  − + − + − + x − + + − Ψ1 (Bm , −Bm ) +A− m σ +Am δ = −gΨ2 (Bm Bm ) + f − Ψ1 (Bβ , −Bβ )ẋ − Aβ σ̇ − Aβ δ̇ where Ψ1 (σ, δ; C1 , C2 ) = C1 cos(σ/2) cos(δ/2) + C2 sin(σ/2) sin(δ/2), Ψ2 (σ, δ; C1 , C2 ) = C1 sin(σ/2) cos(δ/2) + C2 cos(σ/2) sin(δ/2), 4 C1 , C2 ∈ R (7) (8) and 1 (m1 ℓ21 ± m2 ℓ22 ), 4 1 = (ℓ1 f1 ± ℓ2 f2 ). 2 A± m = f± A± β = 1 (β1 ℓ21 ± β2 ℓ22 ), 4 ± Bm = 1 (m1 ℓ1 ± m2 ℓ2 ), 2 Bβ± = 1 (β1 ℓ1 ± β2 ℓ2 ), 2 (9) It is worth noticing that (i) if (3) holds instead of (2), then β0 replaces β after the equals sign in the first equation and each last term of second and third equation (containing ẋ) has to be removed. (ii) In case of identical pendula (m1 = m2 , ℓ1 = ℓ2 ) all the quantities in (9) with superscript − vanish, so that each term in (7) containing them must be cancelled. (iii) Likewise, if also f1 = f2 additional simplifications are evident. (iv) Solving (7) explicitly with respect to the second order derivatives one gets  1h 1 .. + −  Ψ2 (Bm , Bm )(σ̇ 2 + δ̇ 2 )+  Θ x = −kx + g [(m1 + m2 ) sin σ cos δ + (m1 − m2 ) sin δ cos σ] +   2 i  2    1 1  − +    +2Ψ2 (Bm , Bm )σ̇ δ̇ − β0 + 2 (β1 + β2 ) − 2 ((β1 + β2 ) cos σ cos δ + (β2 − β1 ) sin σ sin δ) ẋ,           1 ..  − + + −   Θ σ = −2g mΨ2 (ℓ , −ℓ ) − Ψ1 (−Bm , Bm ) sin δ + 2kΨ1 (ℓ+ , ℓ− )x−   ℓ ℓ 1 2  h i   + − + − 2 2 − +   −Ψ (ℓ , ℓ ) Ψ (B , B )( σ̇ + δ̇ ) + 2Ψ (B , B ) σ̇ δ̇ + 1 2 2 m m    m     m   f1 β1 f2 β2 β2 β1  + −   + + , − + 2βΨ (ℓ , ℓ ) − mΨ − +Θ 1 1  m1 ℓ 1 m2 ℓ 2 m1 ℓ1 m2 ℓ2 m2 ℓ2 m1 ℓ 1 2β − − + + − + −   , Bm )(1 + cos σ cos δ) − Ψ1 (Bm , Bm ) sin σ sin δ] ẋ − Θ(βm σ̇ + βm δ̇), − m [Ψ1 (Bm   ℓ ℓ 1 2          .. 1  + − − +   Ψ (B , −B ) sin δ − 2kΨ1 (ℓ− , ℓ+ )x+ Θ = −2g mΨ (−ℓ , ℓ ) − δ 1 2  m m  ℓ ℓ 1 2  h i   + − − +   +Ψ1 (ℓ− , ℓ+ ) Ψ2 (Bm , Bm )(σ̇ 2 + δ̇ 2 ) + 2Ψ2 (Bm , Bm )σ̇ δ̇ +          β2 f2 β1 β1 β2 f1  − +   − − , + − 2βΨ1 (ℓ , ℓ ) − mΨ1 − +Θ   m1 ℓ 1 m2 ℓ 2 m2 ℓ2 m1 ℓ1  m1 ℓ1 m2 ℓ 2   −  2β   + − − + − +  − m [Ψ1 (Bm , Bm )(1 + cos σ cos δ) − Ψ1 (Bm , −Bm ) sin σ sin δ] ẋ − Θ(βm σ̇ + βm δ̇) ℓ1 ℓ2 (10) 1 where Θ = m − [(m1 + m2 ) + Ψ1 (m1 + m2 , m2 − m1 )] and 2   1 β1 β2 1 ± (ℓ1 ± ℓ2 ), βm = ± . (11) ℓ± = 2ℓ1 ℓ2 2 m1 m2 We notice that the escapement (4) is missing in the first equation governing the beam’s motion. Furthermore, even in the general case of different masses and lenghts the motion of δ (third equation) is more disentangled from the rest of the system than δ: this property will be clearer afterward but actually can be prefigured here, if we imagine to replace in (10) the functions(8) with the second order Taylor polynomial 1 Ψ1 (σ, δ; C1 , C2 ) ≈ C1 − [C1 (σ 2 + δ 2 ) − C2 σδ], 4 2.4 Ψ2 (σ, δ; C1 , C2 ) ≈ 1 (C1 σ + C2 δ). 2 (12) Equilibrium and stability Clearly q = 0 (i. e. x = 0, θ1 = θ2 = 0) is an equilibrium point for system (5) if and only if f1 = f2 = 0 (see (4)) at that position. In that case, since y = 0 if and only if q = 0, i. e. x = 0, σ = 0, δ = 0 is an equilibrium configuration also for system (7) or (10). 5 Whenever F(y) + Φ(y) = 0 (no escapement, no friction), the equilibrium point q = 0 (hence also y = 0) is Lyapunov stable by virtue of the Dirichlet’s criterion, since it is an isolated minimum for 2 P 1 the potential energy V = Kx2 − mj gℓj cos θj . On the other hand, in presence of (2) or (3), the 2 j=1 equilibrium at the same position is asimptotically stable: actually, we can write (2) as D(q)q̇, with ! D= −β −β1 ℓ1 cos θ1 −β2 ℓ2 cos θ2 −β1 ℓ1 cos θ1 −β1 ℓ21 0 −β2 ℓ2 cos θ2 0 −β2 ℓ22 . Since D is a negative–definite matrix, the energy balance d (q̇ · ∇q̇ L − L) = q̇ · Dq̇ < 0 makes the energy a Lyapunov function tending to zero. In a similar way dt one can proceed for the case (3). Nevertheless, if also the escapement (4) is operating, it must be said that this does not connote stability (not even equilibrium) of the system, whichever force F(q) we make use of. 3 Elimination of the escapement. Absence of damping, friction In this paper we are mainly involved in exploring the possibility of in–phase or antiphase synchronization settlements in absence of the escapement (4): from now on, we will investigate the case F(q) = 0 (i. e. no external devices are added to the system). Hence the terms containing f1 , f2 , f pm in (7) or (10) vanish. It must be sais that such as simplification entails a considerable advantage from the mathematical point of view, especially if (4) is a step or discontinuous function, as it occurs in some mentioned models. In our first investigation even the friction forces are temporarily disregarded: the expected fact that synchronization is unattainable will find confirmation. ± Let us now examine the case when also A± β , Bβ (see (9)) are removed: the equations of motion (7) are   d ∇ẏ Le = ∇y Le and, owing to the stability of the configuration y = 0, we are ligitimated now simply dt ! m m1 ℓ1 m2 ℓ2
1 m1 ℓ1 m1 ℓ21 0 and q̇ · Āq̇ − q · V̄ q , where Ā = to replace (1) with the quadratic expansion 2 m2 ℓ2 0 m2 ℓ22 e ẏ), the linearized equations approximating (7) are V̄ = diag (k, m1 ℓ1 g, m2 ℓ2 g). In terms of L(y, .. A1 y +V1 y = 0, A1 = L −1 ĀL −1 = m + Bm − Bm + Bm A+ m A− m − Bm A− m A+ m ! , V1 = L −1 V̄ L −1 = k 0 0 0 0 + Bm g/2 − Bm g/2 − Bm g/2 + Bm g/2 ! or explicitly  .. .. .. − +σ  x +Bm = −kx, +B m δ  m       .. 1 .. .. + − +x − Bm +A+ m σ +Am δ = − g(Bm σ + Bm δ), (13) 2      ..  1 .. ..  + −x − +  Bm +A− m σ +Am δ = − g(Bm σ + Bm δ). 2 Clearly, the same set of equations can be obtained directly from (7), by replacing (8) with the approximations (12) and by neglecting all second order terms. The fundamental frequencies (for both systems in q and y) are found by solving det (λĀ − V̄ ) = 0, leading to   (1 − 2µ)λ3 − λ̄ Y + 2Λ2 (1 − µ + ρ) λ2 + λ̄2 Λ2 (1 + 2γ) λ − λ̄3 Λ2 Y = 0 (14) with µ= m1 + m2 , 2m λ̄ = 2g , ℓ1 + ℓ2 Y = k(ℓ1 + ℓ2 ) , 2mg ℓ1 + ℓ2 Λ= √ , 2 ℓ1 ℓ2 ρ= ℓ1 − ℓ2 m1 − m2 ℓ1 + ℓ2 2m (15) If the two lenghts ℓ1 = ℓ2 = ℓ are the same (but not necessarily the masses m1 and m2 are the same) it results Λ = 1, ρ = 0 and (14) reduces to r g 3 2 2 4 6 (1 − 2µ)λ − ω [Y + 2(1 − µ)]λ + ω (1 + 2Y )λ − ω γ = 0, ω = . (16) ℓ 6 The quantity ω 2 is an evident solution of it: this corresponds to the simplification of the third equation in .. (10) to δ = −gδ/ℓ. The remaining two solutions are (say 1 with − and 2 with +) " 1/2 #  1+Y 4Y 2 2 (17) ω 1± 1− (1 − 2µ) ω1,2 = 2(1 − 2µ) (1 + Y )2 We notice that ω1 6= ω2 , since (1 + Y )2 > 4Y (1 − 2µ), being µ > 0. Moreover, both ω1 and ω2 are different from ω, since µ < 1/2 and from the relation  −1/2 4Y ℓ2 (ω12 − ω 2 )(ω22 − ω 2 ) =B>0 = µℓ 1 − (1 − 2µ) 2g ω12 − ω22 (1 + Y )2 (18) we deduce ω1 < ω, ω2 > ω. Finally, we remark that lim ω1 = ω, µ→0+ p √ lim ω1 = Y ω = k/m, + µ→0 1/2  Y ω, lim ω1 = 1+Y µ→(1/2)− lim ω2 = µ→0+ √ Yω = lim ω2 = ω, µ→0+ lim µ→(1/2)− p k/m, ω2 = +∞ if Y ≥ 1 if 0 < Y < 1 (19) for any Y > 0 (µ ≈ 0 means that the masses of the pendula are negligible with respect to the frame’s one, on the contrary when µ ≈ 1/2 the mass of the frame is negligible). The generalized eigenvectors (λA1 − V1 )v = 0 corresponding to λ = ωj2 , j = 1, 2 and λ̄ are respectively (1, −2(ℓ − g/ωj2 )−1 , 0)T , j = 1, 2 and (0, (1 − m1 /m2 )(1 + m1 /m2 )−1 , 1)T so that the solution y(t) starting from y(0) = (x(0), σ(0), δ(0))T , ẏ(0) = (ẋ(0), σ̇(0), δ̇(0))T is x(t) σ(t) δ(t) x(0) 2ζ2 + C0 ℓ 2 1 + 2 ω1  ẋ(0) 2ζ2 + Ċ0 ℓ 2 #1/2 cos(ω1 t − φ2 ) − 2ζ1 x(0) + C0 ℓ 2 + 1 ω22  2ζ1 ẋ(0) + Ċ0 ℓ 2 #1/2 cos(ω2 t − φ1 ), = B " − B " = 2B ℓ " Y x(0) − ζ1 C0 µℓ 2 1 + 2 ω1  Y ẋ(0) − ζ1 Ċ0 µℓ 2 #1/2 cos(ω1 t − α1 ) − − 2B ℓ " Y x(0) − ζ2 C0 µℓ 2 + 1 ω22  Y ẋ(0) − ζ2 Ċ0 µℓ 2 #1/2 cos(ω2 t − α2 ) − − m1 − m2 m1 + m2 = (20) !1/2 δ̇ 2 (0) δ (0) + cos(ωt − α), ω2 !1/2 δ̇ 2 (0) δ 2 (0) + cos(ωt − α) ω2 2 where B is defined in (18]) and ωj2 δ̇(0) m1 − m2 m1 − m2 δ̇(0), tan α = , C0 = σ(0) + δ(0), Ċ0 = σ̇(0) + ωj2 − ω m1 + m2 m1 + m2 ωδ(0)     ẋ(0) x(0) 1 + Ċ0 / 2ζj + C0 , tan φj = 2 2ζj ωj ℓ ℓ     Y Y 1 ẋ(0) − ζj Ċ0 / x(0) − ζj C0 , tan αj = 2 ωj µℓ µℓ ζj = j = 1, 2 j = 1, 2 We wrote explicitly the solution in order to remark that 1. the motion of δ is independent either of x, σ and of µ, even if m1 6= m2 : as a consequence, none of the initial sets lead to synchronization (i. e. δ(t) → 0), except for the matching of initial data (δ(0) = 0, δ̇(0) = 0): in that case the pendula are exactly in–phase synchronized at any time. 2. The (not null) initial data δ(0) and δ̇(0) produce an effect on the motion of x and σ if and only if the masses m1 and m2 are different. 7 The antiphase synchronization (σ(t) → 0) never occurs, except for the trivial case σ(0) = 0, x(0) = 0, ẋ(0) = 0, σ̇(0) = 0, m1 = m2 (whereas the date for δ are not null, otherwise the system is at rest): for such data σ(t) ≡ 0 (antiphase at each tme) and the frame is at rest. p Y (1 − 2µ) q 3. The solution x(t) is periodic if and only if = , for some rational number q ∈ (0, 1); 1+Y 1 + q2 s if m1 = m2 , then σ(t) is also periodic, otherwise it must be added the condition Y 1 + q2 ∈Q 1 + Y q2 in order to have σ(t) periodic. 4. In order to have some understanding of the amplitudes of oscillation brought along each datum, we represent (18) by means of the parametric functions Bµ (Y ), Y ∈ (0, +∞) with respect to the parameter µ ∈ (0, 1/2): one can easily see that Bµ1 (Y ) < Bµ2 (Y ) if 0 < µ1 < µ2 < 1/2, lim Bµ (Y ) = 0 for any µ ∈ (0, 1/2), Y →+∞ 0 < µ ≤ 1/4 : the global maximum is Bµ (0) = µℓ (strictly r decreasing), ℓ µ 1/4 ≤ µ < 1/2 : the global maximum is Bµ (1 − 4µ) = , 2(1 − 2µ) 2 2Bµ (Y )Y 2 Bµ (Y )ω12 2 Bµ (Y )ω22 (1) (2) , Ψ (Y ) = − and Ψ (Y ) = µ µ µℓ2 ℓ ω12 − ω 2 ℓ ω22 − ω 2 where λj (Y, µ), j = 1, 2, is (17), verify Consequently, the functions Φµ (Y ) = Φµ (0) = 0, lim Φµ (Y ) = 2/ℓ, Φµ1 (Y ) > Φµ2 (Y ) if 0 < µ1 < µ2 < 1/2, 1 1 0 < µ ≤ 1/4 : the global maximum is Ψµ (1/(1 − 4Y )) = p , µ(1 − 2µ) 2ℓ 1/4 ≤ µ < 1/2 : the supremum is 2/ℓ (strictly increasing), (1) (1) (1) (1) Ψµ (0) = 0, Ψµ (1) = 1/2, lim Ψµ (Y ) = 1, Ψµ is strictly increasing for any µ ∈ (0, 1/2), Y →+∞ for any µ ∈ (0, 1/2), Y →+∞ (1) (1) (2) (2) Ψµ1 (Y ) < Ψµ2 (Y ) [resp. >] if 0 < µ1 < µ2 < 1/2 and 0 < Y < 1 [resp. Y > 1], (2) (2) (2) (2) Ψµ (0) = 1, Ψµ (1) = 1/2, lim Ψµ (Y ) = 0, Ψµ is strictly decreasing for any µ ∈ (0, 1/2), Y →+∞ Ψµ1 (Y ) > Ψµ2 (Y ) [resp. <] if 0 < µ1 < µ2 < 1/2 and 0 < Y < 1 [resp. Y > 1]. As for the last point, the quantities Y and µ are considered to be independent of each other by assuming, as an instance, the lenght ℓ and mass of the frame m as fixed and modifying m1 , m2 and k. Finally, we comment the case of different lenghts of the pendula. The circumstance is in whole analogous to the previous case, except for the lesser simplicity of solving (14), in order to achieve expressions similar to (20). With respect to (15), we have Λ = 1 if and only if ℓ1 = ℓ2 , whereas ρ = 0 whether ℓ1 = ℓ2 or m1 = m2 . Hence, Λ > 1 or ρ 6= 0 causes somehow a “disturbance” to the exact solutions we found for ℓ1 = ℓ2 . In regards of that, we only point out that, calling P (λ) the third–degree polynomial of (14), we have P (λ̄) = (1 − Λ2 )(1 − 2µ − Y ) − 2Λ2 ρ. Let us have, for instance, m1 ≤ m2 : assuming ℓ1 < ℓ2 , Y > 1 − 2µ [respectively ℓ1 > ℓ2 , Y < 1 − 2µ ], then P (λ̂) = 0 for λ̂ < λ̄ [resp. >] (see (15)), i. e. the fundamental frequency decreases [resp. increases] by comparison with the case ℓ1 = ℓ2 . Similar remarks can be done about the other two solutions (17), for which it is  2 P (ω1,2 ) = 4 [1 − Y + 2µY (3 + Y )](1 − Λ2 ) + 2ρ[1 − Y 2 + 2Y µ]Λ2 ± (21) o   1/2 ± [(1 + 2µY )(1 − Λ2 ) + 2ρ(1 + Y )Λ2 )] (1 + Y 2 ) − 4Y (1 − 2µ) (22) 4 Elimination of the escapement. Presence of damping, friction Let us start from system (10): defining y = ẋ, u = σ̇, v = δ̇ and setting the system as ẋ = F(x) with x = (x, σ, δ, y, u, v)T , we consider the linear approximation ẋ = (Jx F|x=0 )x, where Jx calculates the Jacobian matrix. Eliminating the terms containing the escapement f1 , f2 and reverting to the explicit expressions of the constant quantities (9), one gets 8                              ẋ = y, σ̇ =  u, δ̇ = v,  1 1 1 −kx + (m1 + m2 )gσ + (m1 − m2 )gδ − β0 y , ẏ = m(1 − 2µ) 2 2 n  g  g 1 (ℓ1 + ℓ2 ) kx − m σ + [m(ℓ1 − ℓ2 ) − 2(m1 ℓ1 − m2 ℓ2 )] δ+ u̇ = 2 2   m(1 − 2µ)ℓ1 ℓ2     β1 β1 β2 β2 + (ℓ1 + ℓ2 )β − m ℓ2 + ℓ1 − − (m1 ℓ1 − m2 ℓ2 ) y − m1  m2  m1 m2   β2 β2 1 β1 1 β1 + − − u− v, 2 m1 m2 2 m1 m2         n    g  g 1   σ − [m(ℓ1 + ℓ2 ) − 2(m1 ℓ1 + m2 ℓ2 )] δ+ (ℓ − ℓ ) −kx + m v̇ =  1 2  2    m(1 − 2µ)ℓ1 ℓ2     2    β1 β1 β2 β2   + −(ℓ1 − ℓ2 )β − m (m1 ℓ1 + m2 ℓ2 ) y − ℓ2 − ℓ1 + −   m1 m2  m1 m2       1 β1 β2 β2 1 β1   − − + u− v 2 m1 m2 2 m1 m2 (23) 2 P 1 (fj,σ σ + fj,δ δ + fj,u u + fj,v v) and mj ℓj 1 ∂fi 1 (f1,σ σ + f1,δ δ + f1,u u + f1,v v) − (f2,σ σ + f2,δ δ + f2,u u + f2,v v) where fi,ζ = , i = 1, 2, m1 ℓ 1 m2 ℓ2 ∂ζ ζ = σ, δ, u, v, are calculated at the equilibrium x = 0, must be added to the fifth and sixth equations, respectively. Remark 4.1 In the presence of escapement (4), the terms j=1 This time, the motion of δ (last equation in (23)) is not independent of the other variables if simply ℓ1 = ℓ2 : actually, this equation is completely disentangled from the rest when also m1 = m2 , β1 = β2 . 4.1 Localization of the eigenvalues The characteristic polynomial associated with the linear system (23) is    β1 β2 β0 6 + λ5 + + (1 − 2µ) (1 − 2µ)λ + m  m1 m2    β1 β2 β0 β1 β2 g m − m1 g m − m2 k + + (1 − 2µ) λ4 + + + + + m2 m1 m2 ℓ2 m  ℓ1  m m  m  m1   g β 2 g m − m2 g k β1 β2 β1 g m − m1 + + + λ3 + + + + m ℓ m ℓ m ℓ m ℓ m m m 1 2 2  1 1 2  2 2 1 g k k g β1 β0 + β2 β2 β0 + β1 g + + λ2 + + + + ℓ ℓ ℓ m m m ℓ m m m 1 1 2  1 22   2 g β g β1 g β2 k g2 k + + + = 0. λ+ ℓ1 ℓ2 m m ℓ2 m1 ℓ1 m2 ℓ1 ℓ2 m (24) The terms with odd exponent of λ are due only to the friction contributions (2): actually, when βj = 0 for each j = 0, 1, 2, the characteristic problem (24) is equivalent to the one formulated in (14). Besides the presence of a Liapunov function which guarantees the asymptotical stability (see Par. 2.1), it is worthwhile to assert also the following Property 4.1 The real part of each root of the polynomial (24) is negative. Proof. It is sufficient to make use of the Routh–Hurwitz criterion (see, for istance, [6]): writing (24) as 6 P g2 k it can be checked, even though calculations last long, that an λn = 0, a6 = 1 − 2µ, . . . , a0 = ℓ1 ℓ2 m n=1 the chain of seven numbers required for the mentioned criterion 1 b2 b1 , a 3 − a 5 , a5   b1 −1   1 a 1 b1 − a 0 a0 b2 a 1 b1 − a 0 b 2 − b1 , a1 − b 2 − b1 , − a0 a5 a3 − a5 a5 a 3 b1 − a 2 b2 b1 b1 a 3 b1 − a 2 b2 a6 , a5 , 9 a0 , where b1 = a4 a5 − a3 a6 , b2 = a2 a5 − a1 a6 , consists of all positive numbers. Since the sequence has no sign change, all the roots of the sixth degree polynomial (24) have negative real parts.  The criterion places the eigenvalues in the left half plane of the complex plane. In order to discriminate the occurrences of real roots, or conjugate pairs of complex roots of (24) one way could be calculate the discriminant which gives additional information on the nature of the roots, real or complex, although calculations for the sixth degree polynomial (24) are quite complex. Our definitive aim is to infer some information about the qualitative behaviour of system (23), by means of locationing as much as possible the portion on the complex plane where the solutions of (24) lie. Remark 4.2 It must be said that the Gershgorin circle theorem used to bound the spectrum is not, in our case, especially powerful: it can be easily seen that the Gershgorin discs where the eigenvalues are confined do not keep them away from the origin of the complex plane: this will be a crucial point in our analysis. 4.1.1 The case of identical pendula Equation (24) will be now discussed in the simpler case of identical pendula: we will assume from now on m1 = m2 = mp , ℓ1 = ℓ2 = ℓ, β1 = β2 = βp (25) (the subscript p is necessary in order not to confuse with the quantities defined in (1) and (2)). Nevertheless, we presume that most of the shown results are still valid in the general case of different pendula. We avoid to write again system (23) in case of assumption (25), since the simplifications are evident. As we touched upon, assumption (25) makes the motion of δ independent of the rest of the system: in fact, .. βp g δ̇ + δ = 0 with the most evident advantage of (25) is the reduction of the equation for δ simply to δ + mp ℓ negative eigenvalues r   p p βp ℓ 1 1 1 βp 2 2 2 2 (26) η ± η − 4 ω if η ≥ 4, − η ± i 4 − η ω if η < 4, η = = − 2 2 ω mp mp g giving the solution  1  !    p p  − ηωt 2 δ̇(0) 1 1  2 2  e 2 η − 4ωt + p η − 4ωt δ(0) cosh sinh    2 2 η 2 − 4ω     1   !  p p − ηωt 1 1 2 δ̇(0) + ηωδ(0) δ(t) = e 2 sin 4 − η 2 ωt + p 4 − η 2 ωt δ(0) cos  2ω  2 2 4 − η     1      − ηωt  1    e 2 δ(0) + δ̇(0) + ηωδ(0) t 2 if η 2 > 4 if η 2 < 4 (27) if η 2 = 4 The characteristic equation (24) can be now written as       βp k g g βp β0 βp β0 2 3 4 λ + λ+ + + (1 − 2µ) + λ + (1 − 2µ)λ + λ2 + ℓ mp m m mp ℓ  m  mp k g β0 + 2βp g k βp λ+ =0 + + m mp m ℓ mℓ (28) where the factorization is related to the uncoupling of δ from the system. We will focus our attention on the factor between square brackets, which gives the eigenvalues related to the motion of x and σ. k ℓ If, in addition to η (see (26)), ω (see (16)), Y = (see (15)), we define mg s β0 ℓ , (29) X= m g then the fourth degree polynomial in square brackets, eq. (28), can be written as (1 − 2µ)λ4 + [X + (1 − 2µ)η]ωλ3 + (ηX + Y + 1)ω 2 λ2 + (ηY + X + 2µη)ω 3 λ + Y ω 4 = 0 −1 We notice that η, X and Y are adimensional quantities (whereas the units of ω are time ). At this point, we employ the Eneström–Kakeya Theorem ([4], [9]), in the following version: 10 (30) Property 4.2 (E–K Theorem) Let pn (λ) = a0 + a1 λ + · · · + an−1 λn−1 + an λn a polynomial with aj > 0 for any j = 0, . . . , n. Then, all the zeros of pn are contained in the annulus of the complex z–plane     a0 a1 a0 a1 an−2 an−1 an−2 an−1 ρm ≤ |z| ≤ ρM , ρm = min , ρM = max (31) , ,..., , , ,..., , a1 a2 an−1 an a1 a2 an−1 an By writing (30) as 4 P an λn = 0 (even though the coefficients are different from the one used in the proof n=1 of Property 4.1), the quantities needed for (31) are Y a0 = ω, a1 X + ηY + 2µη X + ηY + 2µη a1 = ω, a2 ηX + Y + 1 ηX + Y + 1 a2 = ω, a3 X + (1 − 2µ)η X + (1 − 2µ)η a3 = ω. a4 1 − 2µ (32) a2 a1 a3 a0 ≤ and ≤ for any data η, ω, X and Y , hence It is immediate to check that a1 a3 a2 a4     a2 a3 a0 a1 , ρM = max (33) , , ρm = min a1 a2 a3 a4 and the positive quadrant Q = {(X, Y ) | X > 0, Y > 0} is splitted in the four regions     a0 a0 a2 a3 Z1 = (X, Y ) ∈ Q | ρm = , ρM = Z2 = (X, Y ) ∈ Q | ρm = , ρM = a1 a3  a1 a4    a1 a2 a3 a1 Z3 = (X, Y ) ∈ Q | ρm = , ρM = Z4 = (X, Y ) ∈ Q | ρm = , ρM = a2 a3 a2 a4 (34) Sketching a physical reading, we see that a state (X, Y ) ∈ Q on the right part of the quadrant (X ≫ Y ) exhibits a predominance of the damping force on the elastic force, both acting on the beam. On the contrary, in the left upper part of the quadrant (Y ≫ X) the effects are reversed. 4.1.2 Locating the synchronization regions If one opts for the point of view of fixing ω (lenght of pendula) and η (damping of pendula) and let X (friction of the beam) and Y (elastic constant of the spring) vary, one can depict the four regions Zi , i = 1, 2, 3, 4 on the quadrant Q. In finding them, we see that the comparison of (32) leads to the following quadratic conditions in the variables X and Y a0 /a1 < a1 /a2 a0 /a1 < a3 /a4 a1 /a2 < a2 /a3 if and only if if and only if if and only if a2 /a3 < a3 /a4 if and only if X 2 + (η 2 − 1)Y 2 + ηXY + 4µηX + (4µη 2 − 1)Y + 4µ2 η 2 > 0, X 2 + ηXY + X + (1 − 2µ)(η 2 − 1)Y + 2µ(1 − 2µ)η 2 > 0, (η 2 − 1)X 2 + Y 2 + ηXY + ηX + [2 − (1 − 2µ)η 2 ]Y + +1 − 2µ(1 − 2µ)η 2 > 0, X 2 + (1 − 2µ)ηX − (1 − 2µ)Y + (1 − 2µ)[(1 − 2µ)η 2 − 1] > 0 (35) involving the construction of arcs of conics. It is easy to check that for η 2 > 4/3 the first three conditions define three regions on the (X, Y ) positive quadrant delimited by hyperbolae, for η 2 < 4/3 the first and the third conics are real ellipses (η 2 = 4/3: two intersecting lines). The fourth condition refers to a parabola attaining its vertex for some X < 0. The case η ≤ 1, which will be examined deeper, is plotted in Figure 2, where the curves are numbered in the same order as in (35). We make use now of the localization carried out by the E–K Theorem in order to compare the eigenvalues governing δ (see (26)) and those governing σ (solutions of (28)). We will hereafter focus on the case η ≤ 1, which is physically more consistent, asserting that the complementary case can be conceptually treated in the same way. 1 Having in mind (26), we are in the case of conjugate complex eigenvalues for δ with − ηω as real part 2 and ω as modulus. The main question is how the eigenvalues (26) are located with respect to the annulus (31) delimited by the radii (33). We prove the following Proposition 4.1 For η ≤ 1, the two roots (26) cannot lie in the half–plane Re z < −ρM . 11 Y 2 1 Z Z 2 4 1 3 Z 3 Z 4 X Figure 2: The numbered curves are the conics appearing in (35), in the same order. The minimum radius ρm and the maximum radius ρM delimiting the annulus which contains the spectrum are deduced from the inclusion of the status (X, Y) to one of the regions Zi , i = 1, 2, 3, 4, defined in (34). The dotted curves 2 and 3 refer only to the intermediate values in (32) between the minimum and the maximum. 1 Proof. The real part of the two roots (26) is − ηω: they belong to the half–plane Re z < ρM if and 2 ηX + Y + 1 X + (1 − 2µ)η only if (see also (32)) η > 2 in Z1 ∪ Z3 , η > 2 in Z2 ∪ Z4 . However, the two X + (1 − 2µ)η 1 − 2µ 1 1 conditions define empty regions in Q, since they are equivalent to Y + ηX < (1 − 2µ)η 2 − 1 < 0 and 2 2 1 X < − (1 − 2µ)η < 0, respectively (we recall that µ < 1/2, see (15)).  2 Remark 4.3 The eigenvalues (26) belong to the semicircumference |z| = ω, Re z < 0: the more specific question whether they lie in the semicircle |z| ≤ ρM , Re z > 0 can be easily solved, by comparing ω with (33), second couple. It results that (26) are within the semicircle if and only if X > (1 − 2µ)(1 − η) when ρM = a2 /a3 (see (32) and if only if Y > (1 − η)X + (1 − 2µ)η − 1 when ρM = a3 /a4 . Graphically, each of the two regions Z1 ∪ Z3 and Z2 ∪ Z4 (see Figure 4) is splitted by a straight line and the required condition is true only on one side. Our first conclusion is that the system cannot establish a status where the difference δ(t) decays to zero more rapidly than the sum σ(t): thus, the in-phase synchronization onset is inhibited. The result is consistent with the experimental detection, starting from Huygens, on the grounds that the antiphase synchronization is indeed prevailing on the inphase one in this sort of phenomenon. We finally discuss the possibility for the system of establishing a status of antyphase synchronization, still keeping η ≤ 1. Making use once again of the localization (33), we compare the real part of (26) with ρm : whenever ηω < 2ρm , the decay of σ is expected to be faster than the one of δ, so that antyphase synchronization is facilitated. Thus, we are going to check whether (see 32)) (A) η<2 Y X + ηY + 2µη inZ1 ∪ Z2 , (B) η<2 X + ηY + 2µη ηX + Y + 1 inZ3 ∪ Z4 . (36) The two conditions in (36) define two half–planes above the straight lines (2 − η 2 )X + ηY − η(1 − 4µ) = 0 and ηX + (η 2 − 2)Y + 2µη 2 = 0, respectively. As for condition (A) of (36), it can be easily seen that 12 Im z ◊ ◊ ☼ ω ρ M ☼ ω ρ Re z m ☼ ◊ ☼ ◊ Complex plane Figure 3: Thepoints marked with ⋄ are the four solutions of (28), here considered as two couples of conjugate complex roots. The radius ω cannot exceed ρM : two possible values of ω are drawn, the smaller one ω < ρm referring to a state which facilitates the antiphase synchronization. 13 1 2 − η2 then (A) is valid in Z1 ∪ Z2 except for a triangular lower region cut off by the straight 2 4 − η2 line, (see Figure 4) if µ ≤ if 1 2 − η2 1 < µ < then (A) is valid in all the set. 2 24−η 2 By recalling µ = mp /m, we see that, η being equal, the smaller is the mass of the pendula with respect to the mass of the system, the larger is the region which excludes the possibility of antyphase synchronization, i. e. the removed zone. Condition (B) of (36) eliminates all the lower part from Z3 ∪ Z4 and selects the region bounded from below by the straight line ηX + (η 2 − 2)Y + 2µη 2 = 0 and from above by the parabola related to the fourth condition in (35). Moreover if µ ≥ 1 2 − η2 then the selected region is bounded at the left hand side by a segment on the Y –axis, 2 4 − η2 if µ < 1 2 − η2 then the selected region is disjointed from the Y –axis. 2 4 − η2 The different cases are summarized in Figure 4. Y 1 X I Y 1 1 Y X X II III Figure 4: The selection on the quadrant Q is refined by comparing ω with ρm . Picture I refers to the case 1 2 − η2 1 1 1 1 2 − η2 , picture II to < µ < and picture II to ≤ µ < . µ≤ 2 2 24−η 24−η 4 4 2 Let us call A ⊂ Q the subset where conditions (A) and (B) of (36) are both fulfilled. Whenever a state 1 (X, Y ) ∈ A is such that the four roots of (30) are all real, then they exceed in modulus ηω, therefore the 2 decay of δ(t) is lengthier than the decay of σ and the system shows a tendency to an antiphase arrangement. However, if some or all of the roots of (30) are not real, the condition (X, Y ) ∈ A does not guarantee 1 by itself that the real part of such solutions are greater than ηω in modulus. In order to overcome this 2 14 problem, let us explain a method, without delving into the detail: whenever the solutions of (30) are, as an instance, all complex, one can start from the conditions   1 a3 1 a2 2 2 (Re λ) + (37) Re λ + − 2ρm ≥ 0 2 a4 4 a4 where a2 , a3 and a4 are the same as in (32) and Re λ is the real part of any root of (30). The estimation (37) can be proved by setting the complex roots as ξr ± iκr , r = 1, 2, making use of the known relations between roots and coefficients: ξ12 + κ21 + ξ22 + κ22 + 4ξ1 ξ2 = a2 /a4 , 2(ξ1 + ξ2 ) = −a3 /a4 and finally recalling that ξr2 + κ2r ≥ ρ2m , r = 1, 2. In a mixed occurrence of two real roots and two complex conjugate roots one argues in a similar way. Condition (37) can be used in order to separate the real parts away from zero: actually, if we require a2 ≤ 2ρ2m , then Re λ ≤ B, being B the negative solution of (37) where = replaces ≥. Finally, by a4 1 1 demanding − ηω ≥ B, we achieve that the real parts of the solutions of (30) exceed in modulus ηω. 2 2 a2 1 2 ≤ 2ρm and − ηω ≥ B, once they have been expressed in terms of X and Y via (32) The conditions a4 2 and (33), will determine the regions of Q where antiphase synchronization is facilitated. 5 Conclusions In the frame of the study of coupled oscillations of two pendula, we intended to pursue a double objective: 1. to formulate the model in the experimental context as general as possible, 2. to develop the corresponding mathematical problem in a simplified case, drawing the attention to the fact that some classical results about the localization of the spectrum of a matrix can allow us to predict the qualitative behaviour of the system. We selected the parameters X (see (29)) and Y (see (15)) in order to plot on the quadrant Q a certain number of regions in each of which the system develops in a different way. The parameters µ (definrd in (15)) and η (see (26) are considered as constant, but different choices can be made in order to represent the states (as, for instance, fixing the friction of the borad β0 but varying the damping of the pendula βp ). The sections of the complex plane where the spectrum is confined and the regions on Q are not computed via a numerical simulation but they are predicted by the analysis of the spectrum of the linearized system (23). Within specific ranges of the parameters, the tendendy of the system to evolve towards the antiphase synchronization, rather than the inphase one, is predicted by the present analysis. This conclusion is in step with the real development of the phenomenon, starting from the Huygens’ observation in the XV II century of the 180 out of phase swings (see also [5])). Generally speaking, our purpose is to highlight that the method, beyond the specific circumstance which has been exerted to, can be extended to more general situations where the role of certain parameters are exchanged or some restrictions are relaxed. At the same time, the analysis is appropriated, in our mind, to be combined with a simple numerical approach. Both the mentioned points (generalization and matching via computer) are now topics for our current research. More precisely, an in–depth investigation of the problem will concern (i) extending the method to the case of different pendula, by examining equation (24) via the spectrum analysis performed in the simpler case, or by elaborating a formula similar to (21) for evaluating the effects of a variance in the features of the two pendula, (ii) estimating the time of decay of the motion and discard the situations where a very short time from the starting time to the almost rest state would produce not interesting cases, (iii) verifying the analytical outcome by means of simulations via computer, either for calculating the spectrum and for tracing the profiles of δ and σ, (iv) locating the eigenvalues in more restricted regions, by making use of some generalization of the E–K Theorem (as in [8]) which confines the spectrum in specific circular sectors of the complex plane, 15 (v) checking whether the decrease of βi , i = 1, 2, 3 to zero in (23) will lead to the solutions described in Section 3, (vi) adding the effects (4) of an escapement, As for the point (i), the difficulty comes from the non–possibility of factorizing the characteristic equation (24) as in (28), so that the eigenvalues for δ cannot longer be separated from the rest of the spectrum. At this point, the small coefficients (ℓ1 − ℓ2 , . . . ) which join the equation for δ (last equation in (23)) with x and σ will play a significant role. On the other hand, point (vi) renders the analytical problem much complex and deeply different from the present one: as we already remarked, the new formulation requires a non–trivial discussion of equilibrium and stability and of the existence and the regularity of solutions, where the difficulty arises from the typical (but experimentally adequate) discontinuous profile of (4). References [1] Bennett, M. , Schatz, M. F. , Rockwood, H. , Wiesenfeld, K. , Huygens’s cloks, Proc. R. Soc. Lond. A, 458, 563–579 (2002) [2] Czolczynki, K. , Perlikowski, P. , Stefanski, A. , Kapitaniak, T. , Huygens’ odd sympathy experiment revisited, Int J. Bifurcation Chaos 21, 2047 (2011) [3] Dilão, R. , Anti–phase and in–phase synchronization of nonlinear oscillators: The Huygens’s clocks system, Chaos 19, 023118 (2009) [4] Eneström, G. , Remarque sur un théorème relatif aux racines de lequation an xn +an1 xn1 +. . . a1 x+a0 = 0 où tous les coefficientes a sont réels et positifs, Tôhoku Mathematical Journal 18, 34–36 (1920) [5] Fradkov, A. L. , Andrievsky, B. , Synchronization and phase relations in the motion of two–pendulums system, Int. J. of Non–Linear Mechanics, 42, 895–901 (2007) [6] Gantmacher, F. .R. , Applications of the theory of matrices, translated and revised by J. L. Brenner, Interscience Publishers, Inc. , New York (1959) [7] Gelfand, I. M. , Kapranov, M. M. , Zelevinsky, A. V. , Discriminants, resultants and multidimensional determinants, Bulletin (New Series) of the American Mathematical society 37 2, 183–198 (1999) [8] Govil, N. , K. and Rahman, Q. I. , On the Eneström–Kakeya Theorem, Tôhoku Mathematical Journal 20, 126–136 (1968) [9] Kakeya, S. , On the Limits of the Roots of an Algebraic Equation with Positive Coefficients, Tôhoku Mathematical Journal (First Series) 2, 140–142 (1912–13) [10] Kumon, M. , Washizaki, R. , Sato, J. , Kohzawa, R. , Mizumoto, I. , Iwai, Z. , Controlled synchronization of two 1–DOF coupled oscillators, Proceedings of the 15–th Triennal World Congress of IFAC, Barcelon, Spain (2002) [11] Oud, W. , Nijmeijer, H. , Pogromsky, A. , Experimental results on Huygens synchronization, Proceedings of First IFAC Conference on Analysis and Control of Chaotic Systems, Reims, France (2006) 16