Algebraic determination of limit cycles in a family of three-dimensional piecewise linear differential systems

We study a one-parameter family of symmetric piecewise linear differential systems in R which is relevant in control theory. The family, which has some intersection points with the adimensional family of Chua’s circuits, exhibits more than one attractor even when the two matrices defining its dynamics in each zone are stable, in an apparent contradiction with the 3dimensional Kalman’s conjecture. For these systems we characterize algebraically their symmetric periodic orbits and obtain a partial view of the one-parameter unfolding of its triplezero degeneracy. Having at our disposal exact information about periodic orbits of a family of nonlinear systems, which is rather unusual, the analysis allows us to assess the accuracy of the corresponding harmonic balance The first author is partially supported by a MEC/FEDER grant number MTM200803437 and by a CICYT grant number 2009SGR 410. The second and third authors are partially supported by a MEC/FEDER grant number MTM2009-07849. E. Ponce acknowledges the hospitality and support from Centre de Recerca Matemàtica, Bellaterra, Barcelona, Spain, from January to July 2007. ∗Corresponding author Email addresses: jllibre@mat.uab.cat (Jaume Llibre), eponcem@us.es (Enrique Ponce), javieros@us.es (Javier Ros) Preprint submitted to Nonlinear Analysis Series A: Theory and Methods July 10, 2010 This is a preprint of: “Algebraic determination of limit cycles in 3-dimensional piecewise linear differential systems”, Jaume Llibre, Enrique Ponce, Javier Ros, Nonlinear Anal., vol. 74, 6712– 6727, 2011. DOI: [10.1016/j.na.2011.06.051]


Introduction and statement of the main results
In control theory (see [1,2,3,4,5]) an important class of differential systems are the Lur'e systemṡ where A is a n × n constant matrix, b, c ∈ R n , and the input ξ(t) = ϕ(σ(t)) is the feedback of the output σ through a nonlinear continuous function ϕ : R → R. Typically ϕ(0) = 0, so the origin is an equilibrium point. One of the main problems in this theory is to characterize when the origin is a global attractor. Related to this problem there is the Kalman's conjecture [6] which states that if k 1 ≤ ϕ ′ (σ) ≤ k 2 and the linear systemsẋ(t) = (A + kbc T )x(t) have the origin as a global attractor for all k ∈ [k 1 , k 2 ], then the system (1) has also the origin as a global attractor. Kalman's conjecture is true for dimension n ≤ 3 and fails for n > 3, see [7,8].
We remark that these control systems are analytic in R 3 \{Γ − ∪ Γ + } and with Lipschitz constant L = max{ A , B }, we obtain from the standard theorems of ordinary differential equations the existence and uniqueness of a solution for any initial conditions (t 0 , x 0 ) ∈ R × R 3 , and its continuous dependence on initial conditions and parameters. We remark that, in general, the solutions of these piecewise linear control differential systems are C 1 but not C 2 . Note also that if we smooth adequately the characteristic function (2) we can get a monotone smooth characteristic function whose derivative runs through the whole closed interval [0, 1].
We recall that a limit cycle is a periodic orbit isolated in the set of all periodic orbits of the differential system, and that a matrix is Hurwitz or stable if it has all the eigenvalues with negative real part. Using matrices A with −1 as a triple eigenvalue and B with λ < 0 as a triple eigenvalue, Moreno y Suárez [9] showed the existence of one stable limit cycle coexisting with the stable equilibrium point at the origin for a unique fixed value of the parameter λ. Their result was surprising as control practitioners do not expect to find control systems of type (2) having two simultaneous attractors when both matrices A and B are stable. Apparently this also seems to be in contradiction with the 3-dimensional Kalman's conjecture. By computing the eigenvalues of the matrix pencil A + kbc T for k ∈ [0, 1] one can conclude that there is no such contradiction. Later on, in [10] it was shown numerically that there are two intervals of values of λ for which the system studied in [9] has not one but two limit cycles, one stable and the other unstable. One can think that this counter-intuitive coexistence of attractors is related to the triple eigenvalues but this phenomenon is more general as it will be shown in this paper.
Here we will consider system (1) in R 3 satisfying (2) and with e ± = ± (µ 3 + 1, 6µ(1 − µ 2 ), 11µ 2 (1 + µ)) , if µ > 0. From the eigenvalues of A and B µ we deduce that the origin of system (3)- (6) is an attractor if µ < 0, and a repeller if µ > 0, while the two singular points e ± always are attractors for µ > 0. Therefore at µ = 0 we detect a degenerate pitchfork bifurcation of equilibrium points. Note that the situation leads to a triple zero eigenvalue, so that the family of systems under study allows to obtain a specific one-parameter unfolding for such triple-zero bifurcation.
The above eigenvalues have been selected in order to reduce the analytic study of the symmetric periodic orbits to the solving process of a set of polynomial equations. With these values the degree of such polynomials is the lowest possible, but other choices of eigenvalues are possible even they should provide higher degree polynomials and consequently longer computations. In fact such eigenvalue configuration is shown to appear inside the adimensional family of Chua's circuits, see [11] and Section 4.
From (3)-(6) our differential systems are symmetric with respect to the origin of R 3 , i.e. these systems do not change if we change x by −x. In other words if x(t) is a solution of a system (1), then −x(t) is also a solution. When the orbits corresponding to these two solutions coincide, we say that x(t) is a symmetric solution. In what follows we shall study the symmetric periodic solutions of the one-parameter family of differential systems (3)-(6).
Since the differential systemẋ = B µ x is linear in its domain of definition Γ − ∪ S 0 ∪ Γ + , and has no pure imaginary eigenvalues, system (4) has no periodic orbits. Hence it follows that a symmetric periodic orbit intersects the three zones S − , S 0 and S + . In particular a symmetric periodic orbit has at least two points on the plane Γ + . So every symmetric periodic orbit (x(t), y(t), z(t)) can be determined given an initial condition of the form Our main contribution is to provide for the first time in the class of differential systems (3)-(6) an analytical tool, more precisely an algebraic tool, for studying the symmetric periodic orbits. Our algebraic approach is based on the closing equations for symmetric periodic orbits. The use of such closing equations can be traced back to Andronov and coworkers [12] and is also called the Point Transformation method, see for instance [13]. In fact we determine for every value of µ ∈ R all the symmetric periodic orbits having two points in the plane Γ + , as stated in the next result, which constitutes the main contribution of this paper. The possible existence of symmetric periodic orbits with more than two points in Γ + and non-symmetric periodic orbits is to be studied elsewhere.
Theorem 1. Regarding symmetric periodic orbits of system (3)-(6) having exactly 2 points in the plane Γ + the following statements hold.
(a) For µ < µ SN 1 ≈ −17.350 there exist two periodic orbits, collapsing at µ = µ SN 1 to disappear for greater values of the parameter µ. One of these two orbits is stable and the other unstable.
(c) For µ SN 2 < µ < 0 there are two periodic orbits, collapsing at µ = µ SN 2 to disappear for lower values of the parameter µ. One of these two orbits is stable and the other unstable.
From Theorem 1 we see that the stated symmetric periodic orbits are organized in two branches which exist for different ranges of parameter µ.
Moreover, excepting the three values denoted by µ SN 1 , µ SN 2 and µ * , the periodic orbits are hyperbolic when they exist, and so these periodic orbits will persist under sufficiently small perturbations of A, b µ , c and ϕ. At the critical values µ SN 1 and µ SN 2 system (3)-(6) undergoes a saddle-node bifurcation of periodic orbits, that is, two periodic orbit collides to disappear.
Note that the numerical approximations for these critical values can be as precise as desired, since they come from solving algebraic equations. This is rather uncommon in practice when dealing with such global bifurcations.
Theorem 1 is proved in Section 2, where algebraic expressions for determining the mentioned branches of symmetric periodic orbits and their stability are provided.
The determination of non-infinitesimal limit cycles is a very difficult problem which usually is solved numerically or by approximate methods, as the harmonic balance method. The harmonic balance method is a tool for estimating periodic orbits by trying to fit a truncated Fourier series choosing the frequency, amplitude and phase so that the lost of precision comes from the neglected harmonics. Under certain hypotheses, the method can exactly detect Hopf bifurcations, see [14]. Thanks to Theorem 1, we have exact information about periodic orbits so that we can assess the accuracy of harmonic balance predictions.
When the Fourier series takes into account only a single sinusoid, the method is also called the describing function method, see for instance [15,16,17,18]. This method is widely applied in engineering, because it has a beautiful simple interpretation in terms of intersecting graphs. The describing function method can also be used to analyze the dependence of periodic orbits on the parameters of autonomous systems, obtaining a first approach to the bifurcation diagram associated to the periodic orbits of the system: the so-called first harmonic bifurcation diagram, see [18].
Being an approximate method, the predictions of the first harmonic bi-furcation diagram can fail both qualitatively and quantitatively, see [19]. If ϕ : R → R is an odd piecewise linear function with three pieces, then it is shown in [18] for the nonlinear differential control systems of the forṁ all the qualitative behavior of the periodic orbits, this is not the case in dimension 3. In next result we will obtain the first harmonic bifurcation diagram for differential systems (3)-(6) to be compared with the analytic results of Theorem 1. We will show that such diagram is qualitatively but not quantitatively correct. As a consequence, we also shed some light on the non-equivalence in R 3 of Kalman's conditions on the matrix pencil with the global asymptotic stability of the origin, see Remark 1 later.
Proposition 1. Regarding symmetric periodic orbits, the predictions of the describing function method for system (3)-(6) are the following.
(a) For µ < µ ≈ −16.573 two periodic orbits are predicted, collapsing at µ = µ to disappear for greater values of the parameter µ. One of these two orbits is predicted to be stable and the other unstable.
(c) For µ < µ < 0 two periodic orbits are predicted, collapsing at µ = µ to disappear for lower values of the parameter µ. One of these two orbits is predicted to be stable and the other unstable.
(d) For µ ≥ 0 only one periodic orbit is predicted, being stable for µ < 1 and unstable for µ > 1. Proposition 1 is proved in Section 3. In Figs. 2 and 3 we compare the harmonic balance predictions with the exact results provided by Theorem 1 by plotting the amplitude of the different periodic orbits. We have defined A HB = 1 + y 2 HB + z 2 HB , where (1, y HB , z HB ) is the point which satisfiesẋ > 0 at the intersection with the plane x = 1 of the periodic orbit that is predicted by harmonic balance, see [15] and Section 3 for details. Analogously we write

Proof of main results
We start by introducing the algebraic procedure which allows us to convert the problem of determining the symmetric periodic orbits into the problem of solving a set of algebraic equations.

The algebraic procedure
The solution (x(t), y(t), z(t)) of system (5) where and u = e −t . Note that for all positive values of t we will have 0 < u < 1.
Regarding the solution (X(T ), Y (T ), Z(T )) of the system (4)-(6) in the The way for computing the symmetric periodic orbits of system (3)-(6) having two points in the plane Γ + is as follows. Assuming that there exists one of such periodic orbits, let (1, y 0 , z 0 ) ∈ Γ + be the point where this periodic orbit enters into the zone S + ∪ Γ + and let (1, Y 0 , Z 0 ) ∈ Γ + be the point where this periodic orbit exits such zone to enter S 0 . Since this periodic orbit is symmetric it will enter into the zone S − ∪ Γ − through the Equivalently, we can integrate backwards in time the solution from the point We assume in what follows µ = 0, since the case µ = 0 will be treated in a specific way. Then the exponential e −BµT is the matrix If we take Y 0 = y(t) and Z 0 = z(t), then the four unknowns (y 0 , z 0 , t, T ) associated to this symmetric periodic orbit must satisfy the four equations Using the variables u = e −t and v = e −µT instead of t, T, equations (10) become being the e i 's polynomials in their variables.
Note that to every symmetric periodic orbit having two points in the plane Γ + we can associate one solution (y 0 , z 0 , u, v) of equations (11) with It is not obvious (and not necessarily true in general) the converse statement associating every solution (y 0 , z 0 , u, v) of (11) satisfying the above inequality restrictions to one symmetric periodic orbit. However due to the eigenvalues configuration the following result can be shown.
If µ < 0 (µ > 0) and the four conditionsỹ 0 < 0,Ỹ 0 > 0, 0 < u < 1 and v > 1 (0 < v < 1) are fulfilled, then system (3)-(6) has a symmetric periodic orbit passing through the points (1, y 0 , z 0 ) and (1, Proof. Given a solution under the conditions listed, it can be assured that one orbit of system (5) starting at (1, y 0 , z 0 ) enters S + and after a time Similarly we can assure that one orbit of system (4)  Proposition 2 allows to establish an equivalence relationship between solutions of equations (11) and symmetric periodic orbits with two points in Γ + . Once a solution is detected, to study the stability of the corresponding periodic orbit we must locate their characteristic multipliers wit respect to the unit disk of the complex plane. This can be done by using Proposition 3.2 from [20], which assures the similarity between the derivative of the transition map that can be defined from Γ + to Γ − corresponding to the half part of the symmetric periodic orbit and the product of matrix exponentials e BµT e At , when evaluated at the appropriate flight times t and T . Such transition map is the Poincaré half-return map of the periodic orbit. The above product of matrices has always one eigenvalue equal to −1, and if ν 1 and ν 2 are the other two eigenvalues, then ν 2 1 and ν 2 2 are the desired characteristic multipliers of the periodic orbit. In [10] it is proved the following simple algebraic result, useful to avoid explicit eigenvalue computations when we have the trace and determinant of the above product of exponential matrices. Endowed with these tools, the proof of Theorem 1 can be tackled.

Proof of Theorem 1
In order to show Theorem 1, we start by considering the easier case µ = 0, which corresponds with its statement (d) and needs a special treatment.

The case µ = 0
For the solution (x(t), y(t), z(t)) of system (5)-(6) in the region S + ∪ Γ + starting at the point (1, y 0 , z 0 ) when t = 0 we can use (7) and (8). The If we start from (1, y 0 , z 0 ) with y 0 < 0, we see thatẋ = −y 0 > 0 and so we remain for a time t in the right zone until we arrive at (1, Y 0 , Z 0 ), where from (7) with µ = 0 we get (0, Y 0 , Z 0 ) T = e At (0, y 0 , z 0 ) T . Now, using the flow in the central zone and starting from (1, Y 0 , Z 0 ), after a time T we follow the half part of a symmetric periodic orbit if the equations After eliminating Y 0 = −y 0 − T z 0 and Z 0 = −z 0 by using the last two equations, and some standard simplifications we arrive at the algebraic system of equations where we have used (8)  To study the stability of this periodic orbit, we follow the ideas of the end of Section 2.1 and compute the trace and determinant of the matrix e B 0 T e At .
We get We will apply Lemma 1. It is clear from q = −u 6 that |q| < 1. After substituting the values of T and u we see that p = ν 1 +ν 2 = −6(26 √ 3−45) < 0 and 1+q = 30(26 √ 3−45) so that |p| = −p < 1+q clearly holds. Thus ν 1 , ν 2 , and so the multipliers ν 2 1 and ν 2 2 of the periodic orbit, are in the unit circle of the complex plane, so that it is hyperbolic and stable. Thus statement (d) of Theorem 1 is shown.

The case µ = 0
We will analyze the solutions of equations (11) and their relation with the symmetric periodic orbits of system (3) We also note that a dual approach is possible in writing the closing equations. Thus we can define the solution of system (5) in the region S + ∪ Γ + backwards in time starting from the point ( where e −At has the same expression that e At in (8) if we change u by U = 1/u = e t . Then we can take y(t) = y 0 and z(t) = z 0 to write the closing equations in the four unknowns (Y 0 , Z 0 , U, V ) Note that for all t > 0 we will have U > 1, and if µ < 0 (µ > 0) then 0 < V < 1 (V > 1). We get so the following remark. Thus every symmetric periodic orbit must be associated to two different solutions of equations (11).
To simplify a bit equations (11), it is useful to do the translationỹ 0 = y 0 − 6µ andz 0 = z 0 − 11µ 2 in order to reallocate the (y, z)-origin of the plane x = 1. Next, by removing the factor (1 − u)/2 in the equation e 1 = 0, we obtain a simplified first equationẽ 1 = 0, wherẽ The new second equation, after multiplying it by 2µ 2 becomesẽ 2 = 0, wherẽ and to obtain a simplified version of the third equation we define the poly- and a longer expression forz 0 , namely Substituting expressions (14)- (15) in the first and fourth equations and removing the trivial factors −µ/(v − 1) 2 v (u 2 + v 2 ) we obtain two new poly- while the expression of f 2 (u, v) is not shown for sake of brevity. In fact, a more convenient alternative is to take instead of f 2 the polynomial f 1 (u, v) + (1 − u)f 2 (u, v), which after simplifying by v leads to a new polynomialf 2 .
Our next goal will be to characterize the solutions of the polynomial system formed by the two equations From Remark 2, note that if (u, v) are both different from zero and satisfy the polynomial system (18) for a given value of µ, the same is true for the pair (1/u, 1/v). Note also that if we takeỸ 0 = Y 0 − 6µ, we can write for it the same expression that forỹ 0 in (14) by doing the changes u → U and v → V . Then after using the relations u = 1/U, v = 1/V , we get From Proposition 2 and using (14) and (19) we can obtain more quantitative information about the admissible values of (u, v) corresponding to symmetric periodic orbits having exactly 2 points in the plane Γ + .
Lemma 2. Consider the next two sets of inequalities The followings statements hold.  (14) and (19), we arrive, both for positive and negative values of µ, to the two inequalities If µ < 0 then v > 1 and the only possibility to satisfy the first inequality is that v 2 − 4v + 1 = (v − 2) 2 − 3 < 0, and then the second inequality is automatically fulfilled. When µ > 0 we have 0 < v < 1, and now the first inequality holds but the second one requires again v 2 −4v+1 = (v−2) 2 −3 < 0.
Thus statement (a) is true.
For the definition and the basic properties of the resultant of two polynomials see [21] and [22]. Obviously, we can discard all the previous factors so that we will pay attention only to the roots of the equation K(u, v) = 0.
Lemma 3. The number of solutions of the equation K(u, v) = 0 in W , that is, with 0 < u < 1 and 2 − √ 3 < v < 2 + √ 3 is as displayed in Table 1, see also   Table 2.  For the left edge of the rectangle we get so that the unique point that could be a limit point of an interior branch of Range # Solutions # Periodic orbits  namely v 2 and v 5 , see Table 2, satisfying v 2 v 5 = 1. We also obtain where the last two factors have solutions with 0 < u < 1, namely u 0 = 2−  3) verifying (20). After checking the other points we conclude that points P 1 and P 3 are the unique points of the curve K(u, v) = 0 with horizontal Table 2: Numerical values corresponding to marked points in Figure 4. tangent in W , while P 4 is a branching point, since then both K u and K v vanish.
To show the statement, regarding the first two columns of Table 1, it now suffices to take intermediate values ofv in the corresponding range and look for the number of solutions of K(u,v) in (0, 1). This is a standard computation for polynomials in one variable and will be not detailed. with two points in Γ + is as indicated in last column of Table 1.
Moreover the two branches of symmetric periodic solutions with exactly 2 points in the plane Γ + are as indicated in Theorem 1.
Proof. To guarantee that a solution of equations (3)-(6) corresponds to a symmetric periodic orbit with only two points in Γ + , we have to apply Lemma 2 to all the solutions of Lemma 3. We start by checking the sign of µ for such solutions.
From eqs. (16) and (17), we see that the common factors of terms of degree zero in µ are u(v − 1)(v + 1) (v 2 − 4v + 1) . Thus when one follows a branch of solutions shown in Figure 4, that v − 1 is also a common factor of third degree terms in µ of eqs. (16) and (17). In fact, it is not difficult to see that if we denote by µ(u) the value of µ along the branch of solutions from P 3 to P 5 we have lim u→u + 3 µ(u) = −∞, and lim u→u − 5 µ(u) = −∞. Analogously, it can be shown for the branch from P 3 to P 4 that lim u→u − 3 µ(u) = +∞, being µ positive for all the branch from P 3 to P 6 . Then by applying Lemma 2 we discard this branch for counting periodic solutions.
In short, except the branch from P 3 to P 6 , the sign of µ is consistent with Lemma 2. To finish the first statement it remains to show the last inequalities in S 1 and S 2 , depending on the sign of µ.
To check the last inequality in S 1 we start by looking for possible values of (u, v) in the branches P 0 -P 7 or P 3 -P 5 where such inequality is not true any longer, that is Using the above expression to simplify (16) and (17) and removing trivial non-vanishing factors, we obtain the two conditions Eliminating µ between these two equations and between the first of them and equation (21), we get from the polynomial resultants, removing again trivial non-vanishing factors, the conditions Computing now the polynomial resultant of h 1 and h 2 with respect to u, we obtain a polynomial in v whose only roots in the interval (1, and v ≈ 1.267616. Going back with these values to h 1 = 0 and h 2 = 0, we obtain the points (u, v) = (−2, 2) and (u, v) ≈ (0.205425, 1.267616). Clearly, only the last point is in W being in fact on the discarded branch P 3 -P 6 . The conclusion is that, by continuity, the sign of the left hand side of (21) does not change along the branches P 0 -P 7 and P 3 -P 5 .
Since at P 0 the left hand side of (21) is equal to −4, we can assure that the last inequality in S 1 holds for all the points of the branch P 0 -P 7 . We take now an arbitrary point of the branch P 3 -P 5 to check that such inequality is also fulfilled along the branch. We select u = 1/2 and compute the values of Then we conclude from Lemma 2 that all the points of the branches P 0 -P 7 and P 3 -P 5 correspond to periodic orbits of system (3)- (6).
It remains to show that the same is true for the branch P 0 -P 2 , where the sign of µ is positive, and now we must check the last inequality in S 2 .
We start by looking for possible values of (u, v) in the branch where such inequality is not true any longer, that is Using the above expression to simplify (16) and (17) and removing trivial non-vanishing factors, we now obtain the two conditions Eliminating µ between these two equations and between the first of them and equation (22), we get from the polynomial resultants, removing again trivial non-vanishing factors, the conditions Computing now the polynomial resultant ofh 1 andh 2 with respect to u, we obtain a polynomial in v whose only roots in the interval ( Clearly, no point is in W and, by continuity, the sign of the left hand side of (22) does not change along the branch P 0 -P 2 .
Since at P 0 the left hand side of (22) is equal to −4(2 − √ 3) 2 , we can assure that the last inequality in S 2 holds for all the points of the branch P 0 -P 2 . We then conclude from Lemma 2 that all the points of the branch P 0 -P 2 correspond to periodic orbits of system (3) To show the second statement, we first consider the branch from P 3 to P 5 . By considering the common factors of the coefficients of third degree in (16) and (17), we conclude that in both ends the value of µ tends to −∞. By continuity the value of µ(u) reaches a maximum negative value. This value corresponds with µ SN 1 , which can be computed by standard procedures with as much accuracy as desired.
The value of µ in P 1 is positive, being also positive for all points from P 0 to P 2 . Furthermore, reasoning as before, lim u→u − 2 µ(u) = +∞. From P 0 to P 7 , the value of µ is negative and, considering the common factors of the coefficients of zero degree in (16) and (17), it vanishes in both endpoints. By continuity the value of µ reaches a minimum negative value corresponding to µ SN 2 and the statement follows.
From Lemma 3 and Lemma 4, Theorem 1 is shown with the exception on the stability assertions about the periodic orbits. It suffices to compute the trace and determinant of the above product of matrix exponentials and apply Lemma 1. After tedious but standard computations, the stability of the periodic orbits can be obtained and the proof of Theorem 1 is complete.

Computations leading to the first harmonic bifurcation diagram
Before applying the describing function method to our system (3)-(6), we recall that under their assumptions, if σ(t) ≈ σ 1 sin ωt then ϕ(σ(t)) ≈ σ 1 N(σ 1 ) sin ωt, where N(.) is the describing function of the normalized saturation representing the nonlinear gain factor for the first harmonic amplitude, that is N(σ 1 ) = 1 for 0 < σ 1 < 1, and for σ 1 ≥ 1. Thus N(σ 1 ) is decreasing for σ 1 > 1 and tends to zero when and then we arrive at the scalar differential equation in the output σ Thus equation (24) arises by assuming jω to be a root of the characteristic polynomial of the scalar linear substitution differential equation that results by substituting in (25) the nonlinearity ϕ(σ(t)) by N(σ 1 )σ(t), see [3,18] for more details. Note that in our case c T x(t) = x(t), so that for every solution pair (σ 1 , ω) of (24) the describing function method predicts an elliptic periodic orbit given by where we have used ϕ(c T x(t)) ≈ σ 1 N(σ 1 ) sin ωt and (3)-(6).
Typically the determining equation (24) is solved graphically by plotting in the complex plane the curve G(jω) for ω > 0 and the set corresponding to −1/N(σ 1 ), which in our case is the interval (−∞, −1], and looking for possible intersections of both sets. Here we will solve the determining equation (24) in an analytical way as follows. We introduce for any σ 1 ≥ 1 a negative auxiliary variable r = r(σ 1 ) = 1 − 1/N(σ 1 ) ≤ 0, and solve instead of (24) the equivalent equation where the last equality comes from a corollary of the Schur's lemma, see [2]. Now we give a proof of Proposition 1, and next we describe the computations needed to get A HB , the harmonic balance amplitude approximation for the point (1, y 0 , z 0 ), to be compared with the exact value A C .
In short the harmonic balance method predict two periodic orbits for µ < µ and for µ < µ < 0. Consequently the analysis predicts saddle-node bifurcations of periodic orbits at such two negative values of µ. At µ = 0, it is also predicted the disappearance of one periodic orbit (in some kind of heteroclinic connection, since then Ω − = 0). For µ > 0, only one periodic orbit remains, namely the associated to the pair (r − , Ω + ).
The stability of the predicted limit cycles can be deduced in several ways, see [17] or [18]. For the symmetric periodic orbits predicted by harmonic balance in R 3 , the notions of axial stability and radial stability are useful, see [18]. Roughly speaking, the axial stability is associated to the third root ρ of the characteristic polynomial corresponding to the linear substitution problem for every solution (σ 1 , ω) of the determining equation. Such polynomial is in our case which for every solution (σ 1 , ω) of the determining equation has the roots ±jω and ρ. Thus we can deduce easily from the coefficient of s 2 in s 3 + 6s 2 + 11s + 6 − 1 1 − r 6(µ + 1)s 2 + 11(µ 2 − 1)s + 6(µ 3 + 1) that ρ = 6(µ + r)/(1 − r). Since r is negative by definition, sign(ρ) = sign(µ + r) and the axial stability can be guaranteed always for µ < 0.
However for µ > 0 the axial stability can be lost. In fact from (30) and using r − , we can write and this expression becomes positive for µ > 1. Therefore the predicted symmetric periodic orbit for µ > 1 is unstable.
The radial stability can be characterized by computing for every solution of the determining equation the sign of the following expression: and it is required to have positive sign for stable cases. This is equivalent to the so-called Loeb criterion, see [17]. We obtain so for (d/dω) Im G(jω) the where the numerator num(ω) is exactly the left hand side of equation (29) determining the value of ω for the periodic orbit. Denoting such value for ω * and using that num(ω * ) = 0, we see that d dω 6ω num(ω) den(ω) ω=ω * = 6ω * d dω num(ω) den(ω) ω=ω * = 6ω * num ′ (ω * ) den(ω * ) .
The proposition is shown.

Realization in the Chua's circuit
We will see that system (3)-(6) is realized by certain set of parameters in the celebrated Chua's oscillator, see [11] and references therein. The dimensionless equations of the Chua's oscillator arė  See [23] or again [11] for the meaning of different parameters. This system can be written in the form used for system (3)-(6) by the linear change of variables X = x, αY = −(1 + γ)x − y, αZ = (γ 2 − β)x + γy + z, well defined for α = 0.
From the application of Theorem 1 we now deduce that in the first two cases we must expect global asymptotic stability for the origin, while in the third case, as 0 < µ 3 < µ * , the origin is a repeller and there exists a stable limit cycle. In the last case, by using (18)