Unfolding the fold-Hopf bifurcation in piecewise linear continuous differential systems with symmetry

Depósito de Investigación de la Universidad de Sevilla https://idus.us.es/ This is an Accepted Manuscript of an article published by Elsevier in Physica D: Nonlinear Phenomena, Vol. 250, on May 2013, available at: https://doi.org/10.1016/j.physd.2013.01.010 Copyright 2013 Elsevier. En idUS Licencia Creative Commons CC BY-NC-ND Unfolding the fold-Hopf bifurcation in piecewise linear continuous differential systems with symmetry Enrique Ponce∗, Javier Ros, El´ısabet Vela Departamento de Matem´atica Aplicada II, Escuela T´ecnica Superior de Ingenier´ıa, Universidad de Sevilla, Camino de los Descubrimientos, 41092 Sevilla, Spain. Abstract Three-dimensional symmetric piecewise linear differential systems near the conditions corresponding to the fold-Hopf bifurcation for smooth systems are considered. By introducing one small parameter, we study the bifurcation of limit cycles in passing through its critical value, when the three eigenvalues of the linear part at the origin are at the imaginary axis of the complex plane. The simultaneous bifurcation of three limit cycles is proved. Conditions for stability of these limit cycles are provided, and analytical expressions for their period and amplitude are obtained. Finally, we apply the achieved theoretical results to a generalized version of Chua’s circuit, showing that the fold-Hopf bifurcation takes place for a certain range of parameters. Key words: Piecewise linear differential systems, bifurcations, limit cycles ∗Corresponding author Email addresses: [email protected] (Enrique Ponce), [email protected] (Javier Ros), [email protected] (El´ısabet Vela) Preprint submitted to Physica D January 15, 2013 1. Introduction and preliminary results The study of limit cycles of a differential system is, after the analysis of equilibrium solutions, one of the most important problems in the qualitative analysis of dynamical systems. Several tools from bifurcation theory are available in the case of smooth systems in order to guarantee the bifurcation and existence of limit cycles, see [1, 2, 3]. The situation is worse in the case of piecewise smooth systems, see [4], so that new results concerning this class of differential non-smooth systems, which is becoming more and more relevant from the point of view of applications, are needed. By considering for instance the situation associated to the classical Hopf bifurcation for smoth systems, in [5] authors obtained bifurcation results for limit cycles in planar non-smooth systems, see also [6]. For discontinuous cases, we refer to [7, 8, 9, 10]. Later, several results were extended to 3D continuous piecewise linear system, see [11] and [12]. Other results for these systems appeared in [6, 13, 14, 15, 16]. A specific bifurcation in 3D vector fields is the fold-Hopf bifurcation, also called Hopf-zero bifurcation. This bifurcation is characterized by the simultaneous appearance of three eigenvalues at the imaginary axis of the complex plane, see Section 7.4 of [2]. In a recent paper [17], as a first step to study the possible occurrence of an analogous to the fold-Hopf bifurcation in continuous piecewise linear systems with symmetry, authors assumed a linear part controlling the dynamics at the origin near a critical situation with a pure imaginary pair and a single zero eigenvalue. However, strong assumptions on the spectra for the external linear parts involved were required in the quoted paper. Here, we suppress such restrictive assumptions and consider 2 much more general linear parts, maintaining the eigenvalue structure leading to the analogous of a fold-Hopf bifurcation. Thus, we give for the first time information about the general unfolding of such bifurcation in the framework of piecewise linear systems. We consider the following family of piecewise linear differential systems written in the Lur’e form, ˙x =F(x) =ARx+bsat(x),(1) where x= (x, y, z)T∈R3,the saturation function is given by sat(u) =          1 if u > 1, uif |u|⩽1, −1 if u < −1, and the dot represents derivative with respect to the time τ, AR=     t−1 0 m0−1 d0 0     and b=     T−t M−m D−d     ,(2) where t,m,d,T,Mand Dare certain coefficients. System (1)-(2) has the following properties: (a) It is symmetric with respect to the origin, i.e. F(x) = −F(−x). (b) In the region with |x|⩽1 it becomes the homogeneous system ˙x(τ) = ACx=     T−1 0 M0−1 D0 0     x.(3) 3 (c) The coefficients t,m,dand T,M,Dare the linear invariants (trace, sum of principal minors and determinant) of the matrices ARand AC, respectively. Note that AC=AR+beT 1, where e1= (1,0,0)T, and that the considered family of systems is in the so-called generalized Li´enard form, see [18], which is in fact equivalent to the observable canonical form in control theory [19]. Thus, under generic conditions for every system of the form (1) after some change of variables, we can get the matrices in the form given in (2) and (3). Our goal is to study the phenomena associated to the fold-Hopf bifurcation at the origin. To do this, we introduce εas the bifurcation parameter such that the three eigenvalues of matrix ACare −εand ρε ±ωi, where ρ∈Rand ω∈R+are auxiliary fixed parameters. Thus for ε= 0 the three eigenvalues are 0 and ±ωi, which are located on the imaginary axis of the complex plane. Accordingly we choose T(ε) = (2ρ−1)ε, M(ε) = ω2+ρε2(ρ−2), D(ε) = −ε(ρ2ε2+ω2),(4) to be assumed hereafter. We emphasize in the next result an invariant property of systems (1)-(2), whose proof is direct and will be omitted. Proposition 1. System (1)-(2) is invariant under the following transformation (x, y, z, τ, t, m, d, ε)−→ (x, −y, z, −τ, −t, m, −d, −ε). This symmetry property is useful for simplifying the analysis of the family under consideration. Next, we summarize the main assertions about possible equilibrium of the family. 4 Proposition 2. For system (1)-(2) under conditions (4) with ρ6= 0 and ω > 0the following statements hold. (a) If dε < 0the unique equilibrium point is the origin. (b) If dε > 0then the equilibria are the origin and the two points x+ ε=1−D(ε) d, T(ε)−tD(ε) d, M(ε)−mD(ε) dT ,x− ε=−x+ ε. (c) If ε= 0 then all the points of the segment {(x, y, z)T∈R3: (x, y, z)T=µ(1,0, ω2)T,|µ|⩽1} are equilibria for the system. If furthermore d6= 0, the above segment captures all the equilibrium points. Proposition 2 is shown in Section 4. From its statement (c) we see that at ε= 0 system (1)-(2) has a degenerate pitchfork bifurcation, see Figure 1. Note that for dε < 0 the points x± εare vanishing points for the vector field corresponding to |x|>1 but they are out of their corresponding zones. They do not constitute real equilibria, although they still organize the dynamics of external regions. This type of equilibrium is usually called a virtual equilibrium point. In order to study the existence of periodic orbits in system (1)-(2), we start by considering the central zone of system (3) with ε= 0. In this case a direct computation shows that the solution of (3) starting from an arbitrary 5 point (x0, y0, z0)Tfor τ= 0 is given by x(τ) = exp (ACτ)     x0 y0 z0     =1 ω2     ω2cos ωτ −ωsin ωτ 1−cos τ ω3sin ωτ ω2cos ωτ −ωsin ωτ 0 0 ω2          x0 y0 z0     . (5) If |x0|<1 then the corresponding orbit is an arc of the ellipse {(x, y, z)T∈R3:x(τ)−z0 ω22 x0−z0 ω22+y0 ω2 +y2(τ) ω2x0−z0 ω22+y2 0 = 1, z(τ) = z0}, (6) passing through the initial point. For some initial points we get complete ellipses in the region |x| ≤ 1 and then we will have that −ω2< z0< ω2. Therefore, the solutions (5) of system (3) that are contained for all τin the central zone are periodic orbits of system (1)-(2) define a bounded set completely foliated by periodic orbits, which has the shape of two solid cones sharing the elliptic disc ω2x2+y2⩽ω2in the plane z= 0 as their common basis, see Figure 2. The structure of periodic orbits of Figure 2 is independent on the values of t,mand d. For ε6= 0 and small, most of these periodic orbits will disappear but some of them could give rise to limit cycles. In this paper, we analyze the possible appearance of limit cycles from the above set of periodic orbits for ε6= 0 and small. As is well known, a linear system cannot have limit cycles. So that, limit cycles must cross the planar frontiers x= 1 or x=−1. If for some values of parameters there exits a limit cycle living in two zones it will be called bizonal limit cycle. Analogously, we speak of tri-zonal limit cycles when they cross 6 the two separation planes. It should be noticed that, due to the symmetry of the vector field with respect to the origin, non-symmetric limit cycles should come in pairs, that is, when a non-symmetric limit cycle exists then there must be its symmetrical one. This is always the case for bizonal limit cycles. The rest of the paper is organized as follows. Our main results related to the bifurcation and stability of bizonal and tri-zonal limit cycles are stated in Section 2. In Section 3, we show how this non-smooth version of the fold-Hopf bifurcation appears in a generalized version of the Chua’s circuit. Auxiliary results and the proofs of main theorems are provided in Section 4. 2. Statement of main results Our first result concerns the possible bifurcation of symmetrical periodic orbits using the three zones. We introduce a new parameter δ=d−tω2, which characterizes the criticality of the bifurcation, in a similar way to what happens in the cases considered in [11] and [12]. Theorem 3. Let us consider system (1)-(2) under conditions (4) where it is assumed ρ6= 0 and δ=d−tω26= 0 and fixed. For ε= 0 the system (1)-(2) undergoes a tri-zonal limit cycle bifurcation, that is, from the configuration of periodic orbits that exists in the central zone for ε= 0, one limit cycle appears for ρδε > 0and |ε|sufficiently small. It is symmetric with respect to the origin and bifurcates from the ellipse {(x, y, z)T∈R3:ω2x2+y2= ω2,z= 0}. This limit cycle has period P=2π ω+2πρ(ω2−m) ωδ ε−(144π5)1/3tρ5/3[δ2+ω2(m−ω2)2] 5ω5/3δ8/3ε5/3+Oε2, 7 its amplitude in xmeasured as max{x} − min{x}is A= 2 + 3πρω2 2δ2/3 ε2/3− −(π2ρ2ω)2/3[δ(21m+ 2t2−72ω2) + 48(tω2−dm)] 40 ·181/3δ7/3ε4/3+Oε5/3, and it passes through the point P0= (1, y0, z0)T, where y0=−3πρω5 2δ1/3 ε1/3+t 3πρω2 2δ2/3 ε2/3+O(ε), and z0=−πdρω 2δε+3πρω5 2δ1/3 ε4/3+O(ε5/3). Furthermore, the bifurcating limit cycle is stable if and only if t < 0,δ > 0 and d < 0. Theorem 3 is proved in Section 4.1. By using Proposition 1 we could add a new assertion saying that the bifurcating limit cycle is completely unstable (the two characteristic exponents have positive real part) if and only if t > 0, δ < 0 and d > 0. The case δ= 0, that is, d=tω2, when d6= 0 can be shown with the same techniques but requires a special consideration. Theorem 4. Assume for system (1)-(2) under conditions (4) that ρ6= 0, m6=ω2and d=tω26= 0, and fixed. For ε= 0 the system (1)-(2) undergoes a tri-zonal limit cycle bifurcation. For |ε|sufficiently small and ρεt(ω2−m)>0 one limit cycle bifurcates from the ellipse {(x, y, z)T∈R3:ω2x2+y2= ω2,z= 0}. 8 where v1and v2are the voltages across the capacitors C1and C2,i3is the current through the inductance, G=1 R, and the function f(v1) = Gbv1+1 2(Ga−Gb){|v1+E| − |v1−E|} ,(10) models the relevant part of the v−icharacteristics of the nonlinear resistor NL.Eis the saturation voltage of the operational amplifiers in the implementation of the non-linear conductance NL, see Figure 6 (b). This characteristics is responsible for the appearance of multiple equilibria. In fact, after some algebraic manipulation, the existence of non-trivial equilibrium points is equivalent to the existence of solutions for the equation −G 1 + GRN v1=f(v1). The slopes Ga=−1 R6 −1 R3 , Gb=1 R4 −1 R3 , characterize the inner regions where the circuit works, see Figure 7. It must be emphasized that there exist other two more external passive zones represented which are not used at all. In fact, it is usual to represent only the three innermost pieces of the characteristic of Figure 7. As suggested in [20], we may write Chua’s circuit equations (9) in normalized dimensionless form by making the following change of variables: x=v1/E,y=v2/E,z=i3/(EG), and ¯τ=τG/C2. Thus, we get the 15 continuous piecewise linear differential system dx d¯τ=         α[y−bx −(b−a)] , x < −1, α[y−ax],|x| ≤ 1, α[y−bx + (b−a)] , x > 1, dy d¯τ=x−y+z, dz d¯τ=−βy −γz, (11) where a= 1 + Ga/G,b= 1 + Gb/G,α=C2/C1>0, β=C2/(LG2)>0, and γ=RNC2/(LG)<0. In this way, each set of circuit parameters has an equivalent set of five normalized dimensionless parameters {a, b, α, β, γ}. The coefficients of the characteristic polynomials in each linear zone are expressed as follows T=−αa −γ−1, M=α(γ+ 1)a−α+β+γ, D=α[γ−(β+γ)a], t=−αb −γ−1, m=α(γ+ 1)b−α+β+γ, d=α[γ−(β+γ)b], (12) where capital letters correspond to the region with |x|<1. For α6= 0 the so-called observability condition holds, see [18], so that there exists a linear change of variables putting the system in the canonical form (1)-(2). We need to check, apart from the hypotheses of different theorems, the feasibility of conditions (4), that is E1:= (2ρ−1)ε+αa +γ+ 1 = 0, E2:= ω2+ρε2(ρ−2) −α(γ+ 1)a+α−β−γ= 0, E3:= ε(ρ2ε2+ω2) + α[γ−(β+γ)a] = 0. (13) 16 In what follows we will assume b6=a, otherwise the model becomes linear, and we will write α(ε), β(ε) and γ(ε) for certain functions satisfying (13). Taking ε= 0 in the three equations of (13) and looking for positive values of α, we get α0=α(0) = ν−1 a2, β0=β(0) = (1 −a)(1 −a−ν) a2, γ0=γ(0) = 1−a−ν a, (14) where ν2= 1 −a−a2ω2. We see that α0>0 requires ν > 1 and then we have a(1 + aω2)<0, finally arriving at the condition −1/ω2< a < 0. Under this condition it is easy to see that β0>0 and γ0<0 always. From the last three equations of (12) we obtain now t=(b−a)(1 −ν) a2, m=2ν(b−a) + a(2 −a)−b(ν2+ 1) a3, d=(b−a)(1 −ν)(1 −a−ν) a4. Using (14) it is easy to check that the condition det ∂(E1, E2, E3) ∂(α, β, γ)ε=0 =γ0−a(2α0+β0+γ0) + a2α0(1 + γ0)−a3α2 0= = 2ν(1 −ν)/a 6= 0, is fulfilled, so that the Implicit Function Theorem assures the existence of a branch of solutions of (13) starting from the point (α0, β0, γ0), with α(ε)>0, β(ε)>0 and γ(ε)<0 for |ε|sufficiently small. For this set of parameters, it is easy to check for ε= 0 that the non-degeneracy condition δ=d−tω2=−(b−a)ν(1 −ν)2 a46= 0,(15) holds. Then, the appearance of our non-smooth fold-Hopf bifurcation in system (11) for ε= 0 and ρ6= 0 is guaranteed by Theorem 3. The bifurcating 17 tri-zonal limit cycle is always unstable since Sgn(δ) = Sgn(t) = Sgn(d). If we select ρ > 0 all the hypotheses of Theorem 6 are automatically fulfilled and we obtain that two bizonal limit cycles bifurcate and they are stable when a < b. In such case, the three limit cycles appear for ε < 0. From Proposition 2, there also appear two isolated equilibrium points and they turn out to be unstable, as an easy computation shows. 4. Proof of the main results We start by giving the proof of Proposition 2. Proof of Proposition 2. From the definition of system (1)-(2) under conditions (4) in the central zone, it is direct to assure that the unique equilibrium point for ε6= 0 and |x|<1 is the origin. From the third component of equations of system (1)-(2), in the zone x > 1, the x-coordinate of the equilibrium point for ε6= 0 is equal to x+ ε= 1 + ε d(ρ2ε2+ω2). If dε < 0 then x+ ε<1 and the equilibrium point is not a real equilibrium but it is a virtual one. The same is true for the zone x < −1 and x− ε. If dε > 0 then real symmetrical equilibrium points x+ εand x− εare obtained, and the statements (a) and (b) hold. For ε= 0 the computation of the segment of equilibrium points in statement (c) is straightforward. 4.1. Proof of Theorems 3 and 4 In order to prove Theorems 3 and 4 some auxiliary results are necessary. The tri-zonal periodic orbits Γ analyzed in Theorems 3 and 4 are symmetric with respect the origin and they live in the three zones of linearity. 18 Therefore, there exists a point P0= (1, y0, z0)T= (1,p0)Tin Γ located in the plane x= 1, such that when time increases, the orbit will evolve in the zone x > 1, until it transversely crosses the plane x= 1 at the point P1= (1, y1, z1)T= (1,p1)Tafter time τR. Next the orbit goes through the central zone until it hits the plane x=−1 after a time τCat a point P2, which due to the symmetry is −P0. Then Γ goes through the zone x < −1 from the point P2to the point P3=−P1during a time τRand finally, after a time τCit comes back to the point P0, see Figure 8. Note that at the plane x= 1, the contact line y=T(ε) splits the plane in two open half-planes, namely the regions where ˙x > 0 (y < T(ε)) and the one where ˙x < 0 (y < T(ε)). Consequently, we must have y0< T(ε) and y1> T(ε). The symmetry of the system (1)-(2) allows us to determine the periodic orbit Γ by using the half-orbit from P0to P2passing through P1. By the continuity of the flow, it is possible to define in an adequate open set within the plane x= 1 the functions πRand πC, such that πR(p0) = p1 and πC(p1) = −p0so that their composition πH=πC◦πRis the half return map. We denote by τR(p0) the time spent by the orbit from P0to P1, being P0= (1,p0) and P1= (1,p1) points with x= 1 in a neighborhood of P0and P1respectively. Let us call τC(p1) the time employed from P1to P2. See Figure 9. We denote derivatives with respect to the restricted coordinates on the plane x= 1 as Dp(·). Then, if we denote by τ(p0) the total time spent by the orbit to pass from P0to P2, we conclude that Dpτ(p0) = DpτR(p0) + DpτC(p1)DpπR(p0). The eigenvalues of the derivative Dpπat the point p0can be computed by using the next result, which appeared in Proposition 3.2 of [12]. As usual, 19 if a matrix Qis nonsingular, we say that the matrix QMQ−1is similar to matrix M. Proposition 8. Consider a symmetrical tri-zonal periodic orbit of system (1)-(2) starting from (1,p0)and passing through (1,p1)with p0,p1∈R2. Then the product of the two matrices  1DpτC(p1) 0DpπC(p1)  −1DpτR(p0) 0DpπR(p0) = −1Dpτ(p0) 0Dpπ(p0)  is similar to the matrix M= exp (ACτC(p1)) exp (ARτR(p0)) = exp (ACτC) exp (ARτR).(16) Note that the complete return map Π can be obtained by composing π with another half-return map which, due to the symmetry, is identical to π, and so the eigenvalues of Π are the squares of the eigenvalues of π. The next result, which appeared in [13] as Lemma 5, is a straight forward criterion to assure that the two roots of a quadratic equation have modulus less than one. Lemma 9. The two solutions of equation x2−px+q= 0 are inside the unit circle if and only if |q|<1and |p|<1 + q. Using the above lemma, the following remark will be used to determine the stability of period orbits. Remark 10. We note that the matrix Min (16) has one eigenvalue equal to −1. We will denote by λ1and λ2the other two eigenvalues. Therefore, if we 20 take p=λ1+λ2and q=λ1λ2, we have trace(M) = −1+pand det(M) = −q, and thus the characteristic equation of matrix Mis (λ+ 1)(λ2−pλ +q) = 0. Then λ1and λ2are the eigenvalues of the derivative Dpπ, and from Lemma 9 both eigenvalues λ1and λ2are inside the unit circle if and only if |det(M)|<1and |trace(M) + 1|<1−det(M).(17) Finally, as the eigenvalues of the derivative of the complete Poincar´e map Πare λ2 1and λ2 2, if the conditions (17) are fulfilled, we can assure that the corresponding periodic orbit is stable. Before the proof of Theorem 3, we recall a result taken from [5] that will be useful for the proofs. Lemma 11. Let η=ξn(ξ)with nodd, where is a real analytic function in a neighborhood of the origin and such that (0) 6= 0. Then there exists a real analytic function χin a neighborhood of the origin with χ(0) 6= 0 such that ξ=η1/nχ(η1/n). Proof of Theorem 3. In order to determine symmetric tri-zonal limit cycles of system (1)-(2), our strategy will be the following. First, we integrate an orbit that starting from the point P0= (1, y0, z0)Tstays in the zone x > 1, arriving to the point P1= (1, y1, z1)T, after a time τR. If from P0, we integrate the orbit backwards arriving at −P1after a time −τC, then due to symmetry of the problem, we have completed the half of a periodic orbit. This approach allows the elimination of the coordinates y1and z1from the computation of the periodic orbit. 21 The solution of system (1)-(2) in the zone x > 1 passing through the point P0is determined by      x(τ) y(τ) z(τ)     = exp (ARτ)         1 y0 z0     −x+ ε    +x+ ε,(18) where x+ εwas defined in Proposition 2. Analogously, in the central zone, integrating backwards system (1)-(2) from P0one obtains the solution      xb(τ) yb(τ) zb(τ)     = exp (−ACτ)     1 y0 z0     .(19) Obviously expressions (18) and (19) correspond with orbits of system (1)- (2) only if x(τ)⩾1 for all 0 ⩽τ⩽τRand |xb(τ)|⩽1 for all 0 ⩽τ⩽τC. By defining the variable w= (τC, ε, y0, z0, τR) and from equations (18) and (19) we obtain the closing equations b G(w) = 0, where b G1(w) = x(τR)−1, b G2(w) = xb(τC) + 1, b G3(w) = y(τR) + yb(τC), b G4(w) = z(τR) + zb(τC). (20) The system b G(w) = 0constitutes a nonlinear system with four equations and five unknowns, and its solutions with τC>0, τR>0 and y0< T(ε) for |ε|sufficiently small, correspond with periodic orbits symmetrical with respect to the origin. Of course, it has the solution w= (π/ω, 0,0,0,0) which corresponds to the ellipse with z0= 0, y0= 0 and x0= 1 given in (6), namely ω2x2+y2=ω2and z= 0. 22 The closing equations defined by (20) have the spurious solution branch ws= (π/ω, 0, µ, 0,0) for all µ6= 0, that does not correspond to periodic orbits of system (1)-(2) because taking y0=µ6= 0, it is easy to show that max 0<τ<π/ω |xb(τ)|= max 0<τ<π/ω cos(ωτ)−µsin(ωτ) ω>1. This spurious solution branch can be removed by just dividing b G1by a hidden factor τRand using modified closing equations, defined through the functions G1(w) = 1 τRb G1(w), Gi(w) = b Gi(w), i = 2,3,4.(21) The modified closing equations G(w) = 0, have the same solution set of (20) without the spurious solution branch. Choosing τRas bifurcation parameter, we can write the closing equations as G(v, τR) = 0, where w= (v, τR) and v= (τC, ε, y0, z0), and we will parameterize their solutions for vin terms of τR. The corresponding Jacobian matrix is J=DvG(v,0) =         0 2ρ−1−1 0 0πρ ω02 ω2 ω20 0 0 0 0 0 2         ,(22) where v= (π/ω, 0,0,0). Since the determinant of (22) is equal to 2πρω 6= 0, the matrix Jhas full rank and we can apply the Implicit Function Theorem for analytic functions [23], deducing the existence of a function v(τR) = (τC(τR), ε(τR), y0(τR), z0(τR)) in a neighborhood of τR= 0. The corresponding series expansions, were 23 computed through Mathematica [24], as follows τC=π ω−τR+ω2−m 12 τ3 R+ +ω2(15m+ 5t2−9ω2)−6(m2+dt) + mt2 720 τ5 R+O(τ6 R), ε=ω(d−tω2) 12πρ τ3 R+ω[d(9ω2+t2−6m) + t(12mω2−t2ω2−15ω4)] 720πρ τ5 R+O(τ6 R), y0=−ω2 2τR+tω2 12 τ2 R+ω[2(2ρ−1)(tω2−d) + πmρω] 24πρ τ3 R+ +ω2(7mt −6d−t3) 720 τ4 R+O(τ5 R), z0=−dω2 24 τ3 R+ω3(tω2−d) 24πρ τ4 R+O(τ5 R). (23) Defining the parameter δ=d−tω2, and assuming in what follows that δ6= 0, we can invert the series ε(τR) applying Lemma 11, and taking n= 3, ξ=τR,η=ε. We obtain, τR(ε) = 12πρ ωδ 1/3 ε1/3+πρ(δ(9ω2−6m+t2) + 6tω2(m−ω2)) 15ωδ2ε+Oε4/3. (24) The variable τRprovides the time spent by the limit cycle in the zone x > 1. Since this time τRmust be positive for |ε|sufficiently small then the condition ρδε > 0 must hold. Replacing τR(ε) in (23) we obtain the expressions of P= 2 (τC(ε) + τR(ε)), y0(ε) and z0(ε) given in the statement of Theorem 3. 24 (1, y1, z1)Tafter a time τR. If starting from P0, we integrate backwards in the zone x < 1 arriving at P1after a time −τCthen we have completed a periodic orbit. The solution (x(τ), y(τ), z(τ))Tof system (7) in the zone x > 1 passing through the point P0is again determined by (18). In the same way, the backward solution (xb(τ), yb(τ), zb(τ))Tof system (7) starting from P0is given again by (19). By defining the variable w= (τC, ε, y0, z0, τR), we obtain the corresponding closing equations b H(w) = 0, where b H1(w) = x(τR)−1, b H2(w) = xb(τC)−1, b H3(w) = y(τR)−yb(τC), b H4(w) = z(τR)−zb(τC). System b H(ω) = 0constitutes a nonlinear system with four equations and five unknowns, and its solutions with τC>0, τR>0 and y0< T(ε) for ε sufficiently small, correspond with periodic orbits of system (7), whenever x(τ)⩾1 for all 0 ⩽τ⩽τRand xb(τ)⩽1 for 0 ⩽τ⩽τC. The above system has the solutions b w= (2π/ω, 0,0, z, 0) for all z≤ω2, which correspond to ellipses with the same expressions (6) with z0=z, y0= 0 and x0= 1. From all this ellipses, which are tangent to the plane x= 1, we will study the possible bifurcation leading to the existence of limit cycles when ε6= 0 and small. As will be shown, such bifurcation appears from specific ellipse for a certain value z=bz, to be later determined. Apart from the above solutions, for all z≤ω2the closing equations have the spurious solution branch ws= (2π/ω, 0, ν, z, 0) for all z≤ω2and for any ν6= 0, that does not correspond with periodic orbits of system (7) 31 in a similar way to what happens in the proof of Theorem 3. This spurious solutions branch can be removed dividing b H1by τRand defining the modified closing equations, through the functions H1(w) = 1 τRb H1(w), Hi(w) = b Hi(w), i = 2,3,4. Choosing τRas bifurcation parameter, we write the closing equations as H(v, τR) = 0, where v= (τC, ε, y0, z0) and we parameterize their solutions for vin terms of τR. The Jacobian matrix at the point (b v,0) with b v= (2π/ω, 0,0, z), is L=DvH(b v,0) =          0 2ρ−1−1 0 02π ωρ+ 1 ω2z−ρ0 0 ω2−z0 0 0 0−2π ωz0 0          . Since the matrix Lhas not full rank, we can not apply the Implicit Function Theorem and we will use the Lyapunov Schmidt reduction procedure, see [1]. Assuming z6=ω2the rank of matrix Lis equal to three and the image ℑ(L) of the corresponding linear map is generated by the vectors ((1,0,0,0)T,(0,0,1,0)T,0, ρ −ρ+ 1 ω2z, 0, zT), while, ker(L), the kernel of L, is generated by the vector (0,0,0,1)T. We decompose the space R4as a direct sum of the subspace ℑ(L) and its orthogonal complement ℑ(L)⊥, which is generated by the vector (0, ω2z, 0, z(ρ+ 1) −ρω2)T. 32 The orthogonal projection matrix onto ℑ(L)⊥is Qℑ(L)⊥=1 σ2         0 0 0 0 0z2ω40zω2(z(ρ+ 1) −ρω2) 0 0 0 0 0zω2(z(ρ+ 1) −ρω2) 0 (z(ρ+ 1) −ρω2)2         , where σ=pω4z2+ (z(ρ+ 1) −ρω2)2, and the orthogonal projection matrix onto ℑ(L) is obviously Qℑ(L)=I−Qℑ(L)⊥. Note that the vanishing of H(v, τR) implies the vanishing of its two projections onto ℑ(L) and ℑ(L)⊥. The condition Qℑ(L)H(v, τR) = 0leads to         H1 (σ2−z2ω4)H2−zω2(z(ρ+ 1) −ρω2)H4 H3 −zω2(z(ρ+ 1) −ρω2)H2+z2ω4H4         =0,(32) while the condition Qℑ(L)⊥H(v, τR) = 0is         0 z2ω4H2+zω2(z(ρ+ 1) −ρω2)H4 0 zω2(z(ρ+ 1) −ρω2)H2+ (z(ρ+ 1) −ρω2)2H4         =0, and taking into account that the second row of this last expression is linearly dependent on the fourth row, we arrive after reducing non-vanishing factors at the reduced equation g(z, τR) = zω2H2+z(ρ+ 1) −ρω2H4= 0.(33) 33 Analogously, the second row of expressions (32) is linearity dependent on the fourth row, so we can remove the such row and define after removing non-vanishing factor, the vectorial function Φ(v, τR) such that Φ1(v, τR) = H1(v, τR), Φ2(v, τR) = (z(ρ+ 1) −ρω2)H2−zω2H4, Φ3(v, τR) = H3(v, τR), so that system Φ(v, τR) = 0,(34) is equivalent to (32). The decomposition R4= ker(L)⊕ker(L)⊥leads to v= (τC, ε, y0,0) + (0,0,0, z0). We define the new variables τh=τC−2π/ω and zh=z0−z, so that v= (2π/ω, 0,0, z)+(e v, zh) where e v= (τh, ε, y0). In the new variables we can write the system (34) in the form Φ(e v, zh, τR) = 0with Φ(0,0,0) = 0. The Jacobian matrix of Φ evaluated at (e v, zh, τR) = (0,0,0) is equal to De vΦ(0,0,0) =      0 2ρ−1−1 02πσ2 ω20 ω2−z0 0     , and its determinant is different from zero if z6=ω2. Now, the Implicit Function Theorem for analytic functions, see [23], can be applied if bz6=ω2, leading to the existence of the vectorial function e v(zh, τR) = (τh(zh, τR), ε(zh, τR), y0(zh, τR)), such that equation Φ(e v(zh, τR), zh, τR) = 0 is fulfilled. Once the existence of analytic function e v(zh, τR) has been assured we can determine the first terms of its series expansion. 34 These series can be computed assuming undetermined coefficients in the form τC=a10τR+a01zh+a11τRzh+a20τ2 R+a02z2 h+··· ε=b10τR+b01zh+b11τRzh+b20τ2 R+b02z2 h+··· y0=c10τR+c01zh+c11τRzh+c20τ2 R+c02z2 h+··· Substituting these expressions in equation (34) and solving term by term for each order in variables τRand zhwe obtain τC(zh, τR) = 2π ω−τR+1 12τ3 R(ω2−m) + h.o.t., ε(zh, τR) = ω3(ω2−z) [t((ρ+ 1)z−ρω2) + dzω2] 24πσ2τ3 R+ h.o.t., y0(zh, τR) = z−ω2 2τR+(ω2−z)t 12 τ2 R+zh 2τR+ h.o.t. (35) We substitute these series in the reduced equation (33), obtaining g(zh, τR) = 1 12zω2(z−ω2)(−dz −dzρ +tzω2+dρω2)τ3 R+h.o.t.. Defining the parameter δ=d−tω2and assuming in what follows that δ6= 0, the above equation only can have solution for τR6= 0 if we choose as the initial value for ε= 0, the value z=bz, where bz=dρω2 (ρ+ 1)d−tω2=dρω2 dρ +δ, which is well defined and different from ω2, from the hypotheses. Taking now z=bzthe new expression for g(zh, τR) is g(zh, τR) = ω2δ 12 τ3 Rzh−(t+dω2)(dρ +δ)2 12(d2+t2)ρω2τ3 Rz2 h− −dρδω6[5δ+ 2(mt −d)] 240(dρ +δ)2τ5 R+ h.o.t., (36) 35 and for τR6= 0 we can remove the factor τ3 R. Now, since δ6= 0, a new application of the Implicit Function Theorem leads to the existence of a function zh=ζ(τR) such that g(ζ(τR), τR) = 0. The series expansion of such function ζ(τR) has the form ζ(τR) = d1τR+d2τ2 R+O(τ3 R). Introducing this series in (36) and solving for each order in τR, we obtain ζ(τR) = dρω4(3δ+ 2mt −2tω2) 20(dρ +δ)2τ2 R+O(τ3 R). Given that bz=z0−zh, we substitute zh=ζ(τR) in (35) obtaining the following series expansions, τC=2π ω−τR−m−ω2 12 τ3 R+ +t2(m−ω2)−6m2−6δ+ 15mω2−9ω4 720 τ5 R+O(τ6 R), ε=ωδ 24πρτ3 R+ω[6m(δ−tω2)−trδ −9δω2+ 6tω4)] 1440πρ τ5 R+O(τ6 R), y0=−ω2δ 2(dρ +δ)τR+tω2δ 12(dρ +δ)τ2 R+O(τ3 R), z0=bz+dρω4[3δ+ 2t(m−ω2)] 20(dρ +δ)2τ2 R+O(τ3 R). (37) Since δ6= 0, we can invert the series ε(τR) applying Lemma 11, by taking 36 n= 3, ξ=τRand η=ε. We obtain, τR(ε) = 2(3π)1/3ρ ωδ1/3ε1/3+ +2πρ [δ(t2−6m+ 9ω2) + 6tω2(m−1)] 15ωδ2ε+O(ε4/3). (38) Since the variable τRprovides the time spent by the limit cycle in zone x > 1, this time τRmust be positive and the condition ρδε > 0 must hold. By replacing τR(ε) in expressions of (37) we obtain the expressions τC(ε), y0(ε) and z0(ε). The peak-to-peak amplitude of the limit cycle is measured as before, but now we cannot exploit the symmetry. Let us denote by τRmax the time for which x(τ) reaches its maximum value x(τRmax) in the zone with x > 1, and by τRmin the time for which x(τ) reaches its minimum value x(τRmin) in the zone with x < 1, so that, the amplitude is equal to x(τRmax)−x(τRmin). Thus, equalities (25) and (26) apply to this case, and we can also compute a series expansion in the variable τRfor τRmax. We get τRmax =τR 2+t 24τ2 R+O(τ3 R).(39) Using (37) and (39) in (26) we obtain the expansion x(τRmax) = 1 + O(τ3 R). Analogously, the time τRmin has a series expansion in the variable τR. Since the variable xreaches its minimum value for τR=τRmin, the following equality must hold ˙x=T(ε)x(τRmin)−y(τRmin) = 0,(40) 37 where      x(τRmin) y(τRmin) z(τRmin)     = exp (−τRminAC)     1 y0(τR) z0(τR)     .(41) From equalities (37), (40) and (41), the series expansion of τRmin is, τRmin =π ω−τR 2+t 12τ2 R+O(τ3 R).(42) Using (37) and (42) in (41) we obtain the expansion, x(τRmin) = dρ −δ dρ +δ−dρω2[2(d−mt)−5δ] 10(dρ +δ)2τ2 R+O(τ3 R). Substituting τRby τR(ε) from (38) and computing A=x(τRmax)−x(τRmin), we arrive at the expression of Theorem 5 for the amplitude. Using that the period of the orbit is P=τC+τRand (37), we easily get its series. In order to determine the stability of the periodic orbit, we follow Remark 13, and we compute the determinant and trace of matrix Musing (30) to obtain trace(M) = 3 + tτR+t2 2τ2 R+δ(2ρ−1) + 2t3ρ 12ρτ3 R+ +2πd [ρ(2tω2−δ)−δ] + dω3(1 −2ρ) 24πρω2τ4 R+ +t(πtρ(t2−4ω2) + (2ρ−1)ω3)) 24πρ τ4 R+O(τ5 R), (43) 38 and det(M) = 1 + tτR+t2 2τ2 R+δ(2ρ−1) + 2t3ρ 12ρτ3 R+ +πt(t3ρ+ 4ρδ −2δ) + δω(1 −2ρ) 24πρ τ4 R+O(τ5 R). (44) We distinguish first the case t6= 0. For εsufficiently small, we have τR(ε) also small and positive, so that |det(M)|<1 if and only if t < 0. Using (43) and (44) in the additional required condition |trace(M)−1|<1 + det(M), we get the inequality δ(dρ +δ) ρ>0.(45) The hypothesis bz < ω2implies dρ dρ +δ<1. Now, we distinguish two cases. Assume first that dρ +δ < 0. Then dρ > dρ +δand so δ < 0, and then (45) is equivalent to ρ > 0. If now we consider dρ +δ > 0 then we have dρ < dρ +δand so δ > 0. Again we get as before the equivalence of (45) and the condition ρ > 0. Now, if we consider t= 0, then δ=dand from bz < ω2we obtain ρ ρ+ 1 <1, which is equivalent to ρ > −1. We analyzed first the inequality |trace(M)−1|<1 + det(M) what it is equivalent to ρ(ρ+ 1) >0. Since ρ > −1 we need ρ > 0 for stability. If we consider now the condition for stability |det(M)|<1, we will require dρ(2ρ−1) <0, which gives rise to the inequality d(2ρ−1) <0 provided that ρ > 0. Proof of Theorem 6. Since 0 <bz < ω2, the proof of Theorem 6 follows directly from the comments before the statement of Theorem 5. 39 5. Acknowledgements Authors are partially supported by the Ministerio de Ciencia y Tecnolog´ıa, Plan Nacional I+D+I, in the frame of projects MTM2009-07849, MTM2010-20907, MTM2012-31821, and by the Consejer´ıa de Educaci´on y Ciencia de la Junta de Andaluc´ıa under grants TIC-0130 and P08-FQM- 03770. References [1] M. Golubitsky, D. G. Schaeffer, Singularities and groups in bifurcation theory., Vol. 69 of Applied Mathematical Sciences, Springer-Verlag, New York, 1988. [2] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields., Vol. 42 of Applied Mathematical Sciences, Springer-Verlag, New York, 1983. [3] Y. A. Kuznetsov, Elements of Applied Bifurcation Theory., 3rd Edition, Vol. 112 of Applied Mathematical Sciences, Springer-Verlag, New York, 2004. [4] M. Di Bernardo, C. Budd, A. R. Champneys, P. Kowalczyk, Piecewisesmooth dynamical systems: theory and applications, Vol. 163 of Appl. Math. Sci., Springer-Verlag, 2007. [5] E. Freire, E. Ponce, J. Ros, Limit cycle bifurcation from center in symmetric piecewise-linear systems, International Journal of Bifurcation and Chaos 9 (1999) 895–907. 40 Figure 4: Limit cycles of system (1)-(2) under conditions (4) determined by Theorem 3 and Theorem 6 that exist for ρ= 1, ε=−0.01, ω= 1, t=−2, m= 2 and d=−3. The tri-zonal limit cycle which appears in the center of the figure is unstable, and the bizonal symmetrical ones are stable. On the double cone, there appear in thin line the ellipses where limit cycles have bifurcated from. The equilibria of external zones are stable. 47 0.6 0.8 1.0 -0.2 0.0 0.2 0.8 0.9 1.0 1.1 x y z Figure 5: The periodic orbits in the vicinity of the upper equilibrium point for ρ= 1, ε=−0.01, ω= 1, t=−2, m= 2 and d=−3, as in Figure 4. We show a new unstable periodic orbit, numerically detected which is between the stable equilibrium and the stable periodic orbit predicted by Theorem 6. 48 − V− V + V+ V R6 R4 R5 R3 R1 R2 (b) R R0 L C2 C1 R N i (a) 000000 000000 000000 000000 000000 000000 000000 000000 111111 111111 111111 111111 111111 111111 111111 111111 000 000 000 000 000 111 111 111 111 111 00 00 00 00 11 11 11 11 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 000 000 000 111 111 111 000 000 000 000 000 000 111 111 111 111 111 111 0000 0000 0000 0000 1111 1111 1111 1111 000 000 000 000 000 111 111 111 111 111 0000 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 00 00 00 00 00 00 00 00 11 11 11 11 11 11 11 11 000 000 000 000 000 000 111 111 111 111 111 111 000 000 000 000 000 000 111 111 111 111 111 111 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 RN L Ri NL C1 C2 R1 R2 R3 R4 R5 R6 Figure 6: (a) Generalized version of Chua’s circuit, including an extra negative resistance device RNapart from the nonlinear resistance device NL. (b) Detailed circuit of the NL device. 49 v i Ga E −E Gb −G/(1+GR0) −G 1 + GRN −E Ga E Gb Figure 7: Current-voltage characteristic of the nonlinear resistance NL and the graph of the straight line i=−GV 1 + GRN . 50 τC τE P0 x=−1 x=1 1 P P3 2 P P2 τC P1 τR P0 P3 x=−1x= 1 Figure 8: Tri-zonal periodic orbit of system (1)-(2) and the relevant points necessary to build the closing equations. 51 τC τE P0 x=−1 x=1 1 P P3 2 P P2P1 P0 x=−1x= 1 Figure 9: Half-return map of system (1)-(2) from section x= 1 to section x=−1. 52 τC τE P0 x=1 1 P τC P1 P0 τR x= 1 Figure 10: Bizonal periodic orbit of system (7) and the relevant points necessary to build the closing equations. 53