The Focus-Center-Limit Cycle Bifurcation in Symmetric 3D Piecewise Linear Systems

SIAM J. APPL. MATH.c 2005 Society for Industrial and Applied Mathematics Vol. 65, No. 6, pp. 1933–1951 THE FOCUS-CENTER-LIMIT CYCLE BIFURCATION IN SYMMETRIC 3D PIECEWISE LINEAR SYSTEMS∗ EMILIO FREIRE†, ENRIQUE PONCE†,AND JAVIER ROS† Abstract. The birth of limit cycles in 3D (three-dimensional) piecewise linear systems for the relevant case of symmetrical oscillators is considered. A technique already used by the authors in planar systems is extended to cope with 3D systems, where a greater complexity is involved. Under some given nondegeneracy conditions, the corresponding theorem characterizing the bifurcation is stated. In terms of the deviation from the critical value of the bifurcation parameter, expressions in the form of power series for the period, amplitude, and the characteristic multipliers of the bifurcating limit cycle are also obtained. The results are applied to accurately predict the birth of symmetrical periodic oscillations in a 3D electronic circuit genealogically related to the classical Van der Pol oscillator. Key words. piecewise linear systems, bifurcation theory, limit cycles AMS subject classifications. 37G15, 34C15 DOI. 10.1137/040606107 1. Introduction and main results. Piecewise linear modeling of nonlinear dynamical systems is especially successful in some engineering problems, such as the analysis and design of electronic oscillators or control systems (see, e.g., [CFPT02]). However, in the framework of piecewise linear systems, there are no general bifurcation results explaining the appearance or disappearance of self-sustained oscillations, as is the case for the Hopf bifurcation theorem in the context of differentiable systems. Thus, the authors gave in [FPR99] a complete characterization of the focus-centerlimit cycle bifurcation for symmetric planar piecewise linear systems. Now we show how the corresponding result can be extended to the 3D case. We consider a common situation in applications, namely, dynamical systems defined by piecewise continuous vector fields with three linear zones and two parallel frontiers. Furthermore, it is assumed that such systems show symmetry with respect to the origin; that is, if we put them in the form dx/dτ =f(x) with x∈R3, they satisfy f(−x)=−f(x). In particular, f(0) = 0, and so the origin is an equilibrium point for all values of the parameters. By means of a linear change of variables, it is always possible to suppose that the frontiers are the planes Σ1={x∈R3:x1=1} and Σ−1={x∈R3:x1=−1}. We denote by L(left), C(central), and R(right) the regions of R3at which x1<−1, |x1|≤1, and x1>1, respectively, hold. To be more precise, we consider systems expressed as follows: ˙x =⎧ ⎨ ⎩ ALx−bif x1<−1, ACxif |x1|≤1, ALx+bif x1>1, (1.1) ∗Received by the editors April 1, 2004; accepted for publication (in revised form) February 2, 2005; published electronically August 3, 2005. This research was partially supported by grants DPI2000-1218-C04-04, BFM2001-2668, and BFM2003-00336 of the Spanish Ministry of Science and Technology. The authors were also supported by Junta de Andaluc´ıa, as members of the research group TIC-130, in the budget corresponding to 2003. http://www.siam.org/journals/siap/65-6/60610.html †Departamento Matem´atica Aplicada II, Universidad de Sevilla, Escuela Superior de Ingenierosm, Camino de los Descubrimientos s/n, 41092, SEVILLA, Spain ([email protected], [email protected], [email protected]). 1933 Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1934 EMILIO FREIRE, ENRIQUE PONCE, AND JAVIER ROS T>TT<T cT=T cc Fig. 1.The focus-center-limit cycle bifurcation in the case D>0,γ>0. The focal plane and the complementary one-dimensional invariant manifold at the origin are shown, along with the two parallel planes which separate the three linear regions. In the situation sketched, as deduced from Theorem 1.1, the bifurcating limit cycle is of saddle type. where we have taken advantage of the continuity and symmetry of the vector field involved; in particular, the matrices ALand ACdiffer only in their first columns. From Proposition 16 of [CFPT02], under the generic condition of observability, every system (1.1) can be written in the generalized Li´enard form d dτ ⎡ ⎣ x1 x2 x3⎤ ⎦=⎡ ⎣ t−10 m0−1 d00 ⎤ ⎦⎡ ⎣ x1 x2 x3⎤ ⎦+⎡ ⎣ T−t M−m D−d⎤ ⎦sat(x1),(1.2) where sat(x1)isthenormalized saturation sat(x1)=⎧ ⎨ ⎩ −1,x 1≤−1, x1,|x1|<1, 1,x 1≥1, so that, regarding system (1.1), we have AL=⎡ ⎣ t−10 m0−1 d00 ⎤ ⎦,A C=⎡ ⎣ T−10 M0−1 D00 ⎤ ⎦,b=⎡ ⎣ T−t M−m D−d⎤ ⎦. Note that system (1.2) is a particular instance of the more general Lur’e form dx dτ =Ax+bsat(cTx) for the case A=ALand c=e1, where e1stands for the first vector of the canonical basis. Clearly, the parameters t,m,dand T,M,Dstand for the trace, the sum of principal minors of order two, and the determinant of each matrix, and they completely determine the dynamics of the system. Choosing Tas the bifurcation parameter, for the critical value Tc=D/M with M>0, system (1.2) has a linear center in the zone C(see Figure 1); that is, the matrix AChas a pair of pure imaginary eigenvalues. We want to analyze whether a limit cycle bifurcates from this configuration as the bifurcation parameter Tvaries. Note the similarities with the classical Hopf bifurcation scenario. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 3D FOCUS-CENTER-LIMIT CYCLE BIFURCATION 1935 It will be useful, in order to know the stability of such a limit cycle, to estimate the characteristic multipliers of the limit cycle, that is, the eigenvalues of the derivative of a Poincar´e return map defined in an adequate section of the phase space. We will denote the logarithms of these characteristic multipliers by μrand μa, from radial and axial, respectively. Our main result is the following. Theorem 1.1. Let us consider system (1.2) with M>0,Tc=D/M, and γ=DM −Dm +dM −tM2=0.ForT=Tcthe system undergoes a focus-centerlimit cycle bifurcation; that is, from the lineal center configuration in the central zone, which exists for T=Tc, one limit cycle appears for γ(T−Tc)>0and T−Tc sufficiently small. The amplitude “a” (measured as the maximum of |x1|), the period Per, and the logarithms of characteristic multipliers μrand μaof the periodic orbit are analytic functions at 0, in the variable (T−Tc)1/3; namely, a=1+(6π)2/3M4/3 8γ2/3(T−Tc)2/3+(6π4)1/3a4 960M1/3γ7/3(T−Tc)4/3+O(T−Tc)5/3, Per =2π √M+π(M−m)√M γ(T−Tc)−62/3π5/3M5/6P5 20γ8/3(T−Tc)5/3+O(T−Tc)2, μr=−(48π)1/3M7/6γ2/3 D2+M3(T−Tc)1/3+O(T−Tc)2/3, μa=2πD M3/2+(48π)1/3 M5/6Mt−D γ1/3+M2γ2/3 D2+M3(T−Tc)1/3+O(T−Tc)2/3, where a4=−120tM5+120D+2t3+21mt +72dM4 +−93m+27t2D+27m−2t2dM3+2t2m+25dt −27m2DM2 +25D3+23(mt −d)D2M−25mD3, P5=[M(M−m)2+(Mt−d)2](Mt−D). In particular, if γ>0and D<0, then the limit cycle bifurcates for T>T cand is orbitally asymptotically stable. This theorem describes a codimension-one bifurcation, similar to the Hopf bifurcation of differentiable dynamics (see [CH82]), but some differences should be noted. In particular, the expressions characterizing the bifurcation are in terms of the parameter to the one third power instead of the one half power, and, more important, the limit cycle amplitude’s leading order is O(1). Thus, the stability change of the origin is accompanied by the abrupt appearance of a limit cycle of significant size. This comment also applies to the planar case, as appeared in [Kr87] and [FPR99]. When the coefficient γis not equal to zero, it allows a complete characterization of the bifurcation criticality. Its role is analogous to the coefficient of the cubic term in the Poincar´e–Andronov–Hopf normal form. When γ= 0, the bifurcation is of higher codimension, requiring a specific treatment that will appear elsewhere. We want to remark that it is possible, with the same techniques, to obtain similar bifurcation results for the asymmetric case of single-sided saturation. Thus, the proposed methodology is able to cope with a wider class of piecewise linear systems. The rest of the paper is structured as follows. In section 2, we show how the above result can be useful for accurately predicting the birth of symmetrical periodic Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1936 EMILIO FREIRE, ENRIQUE PONCE, AND JAVIER ROS oscillations in a tridimensional electronic circuit, which can be built by taking a Van der Pol oscillator as starting point. The proof of Theorem 1.1 is given in section 3. 2. Predicting the onset of symmetrical periodic oscillations in a 3D electronic circuit. In this section, we consider the electronic circuit of Figure 2(a), genealogically related with the classical Van der Pol oscillator, in order to show the applicability of our results. Regarding this circuit, the nonlinear conductance NL is its active element, implemented by means of an operational amplifier with the feedback structure of Figure 2(b), and the current-voltage characteristic is shown in Figure 2(c). Note that we are dealing with a nonlinearity characteristic qualitatively similar to the cubic one appearing in the classical Rayleigh–Van der Pol oscillator. In fact, if we eliminate the capacitor C2and make R=R0= 0, then the resulting planar circuit could be thought of as a modern electronic realization of such classical oscillators; see [Kr87] and [FPR99]. Thus, the 3D circuit of Figure 2 can be built by adding the capacitance C2to a bidimensional oscillator circuit. In the context of chaotic circuits, such topology was originally proposed in [SYM81], and was studied afterwards in [FGA84] and [FRGP93] in the case R0= 0 and assuming a nonlinear positive conductance for the resistor R. With slight modifications, this circuit has been extensively studied in the last two decades; see [GK92] or [HBCJM91]. Taking R0= 0 and substituting the nonlinear element by the so-called Chua diode, many papers have also been written; see [CWHZ93], [Ma93], and references therein. Anyway, the onset of symmetrical periodic oscillations was never accurately predicted, since in most cases the circuit was analyzed by taking polynomial approximations. Thus, the rapid bifurcation for the limit cycle observed in practice was never justified. It should be remarked that the characteristic of Chua’s diode is qualitatively similar to the one presented in Figure 2(c) but the zone of negative slope is made up by three pieces with two different slopes. For that, at least two subcircuits with operational amplifiers like those shown in Figure 2(b) are needed. Thus, the Chua circuit characteristic has five linear segments instead of only three, as in our case. However, in modeling Chua’s circuit, usually only the three innermost pieces are represented, since the two outermost pieces of positive slope are not physically used; see [Ke93]. As stated in [Kr87] and [FPR99], there exists an excellent agreement between the actual response of the nonlinear device NL in the circuit and its symmetric piecewise C2 R L R0 NL R3 R1 R2 i v (b) (c) C1 (a) iLi Fig. 2. (a) The 3D electronic circuit. (b) Implementation of the nonlinear conductance NL. (c) Piecewise linear current-voltage characteristic of NL. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 3D FOCUS-CENTER-LIMIT CYCLE BIFURCATION 1937 linear mathematical model. Therefore, we are led to consider the piecewise linear dynamical system C1 dv1 dτ =v2−v1 R−i(v1), C2 dv2 dτ =v1−v2 R−iL,(2.1) LdiL dτ =v2−R0iL, where v1and v2are the voltages across the capacitors C1and C2, respectively, while iLis the current through the inductance. The nonlinear current-voltage characteristic is i(v1)=v1−f(v1) R1 with f(v1)=Esign (v1),|v1|> E/σ, σv1,|v1|≤E/σ, where σ=1+R2 R3 is the gain of the operational amplifier configured (using feedback) as a noninverting amplifier and Eis its saturation voltage. With the following linear change of variables and time rescaling, v1=E σy1,v 2=E σy2,i L=E σC2 Ly3,τ=RC1¯τ,(2.2) and defining the following five nonnegative dimensionless parameters, r=R R1 ,c=C1 C2 ,μ=(σ−1) R R1 =RR2 R1R3 ,ρ=R2C2 1 LC2 ,κ=RR0C1 L,(2.3) we can express system (2.1) as follows: d d¯τ⎡ ⎣ y1 y2 y3⎤ ⎦=⎡ ⎣−r−11 0 c−c−√ρ 0√ρ−κ⎤ ⎦⎡ ⎣ y1 y2 y3⎤ ⎦+⎡ ⎣ μ+r 0 0⎤ ⎦sat (y1).(2.4) For the subsequent analysis, we will choose μand ρas the main bifurcation parameters. In practice, to detect the bifurcation in a experimental way, it is better to tune the parameter μby means of a variable resistor R2, which is equivalent to varying the gain σ. The observability matrix for system (2.1) is O=⎡ ⎣ eT 1 eT 1A eT 1A2⎤ ⎦=⎡ ⎣ 100 −r−110 (r+1) 2+c−c−r−1−√ρ⎤ ⎦, which has full rank for all the values of components of the circuit. From Proposition 16 of [CFPT02], system (2.1) can be expressed in Li´enard’s generalized form (1.2) with the following values: T=μ−c−κ−1,t=−r−c−κ−1<0, M=(c+1)κ−(c+κ)μ+ρ, m =(c+1)κ+(c+κ)r+ρ>0, D=(cκ +ρ)μ−ρ, d =−(cκ +ρ)r−ρ<0. (2.5) Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1938 EMILIO FREIRE, ENRIQUE PONCE, AND JAVIER ROS M = 0 D = 0 MT = D –0.05 0 0.05 0.1 0.15 0.2 0.25 ρ 0.2 0.4 0.6 0.8 1 μ Fig. 3.The parabolic arc (thick line) in the plane (μ, ρ)corresponding to the bifurcation locus of Proposition 2.1 for c=0.2and κ=0.05. The horizontal line indicates the path followed as μ varies for a fixed value of ρ. The dashed line represents points with D=0, so that above it we have D<0. At the dotted straight line we have M=0, and above this line we have M>0. The vertical line corresponds to μ=μ∗. Note that these coefficients are the linear invariants of the two matrices involved, so that their computation is straightforward, and that it is not necessary to explicitly compute the linear change of variables required to get the Li´enard form for applying Theorem 1.1. The equation MT −D= 0 leads to (c+κ)μ2−(c+κ)2+c+2κμ+ρ(c+κ)+κ(c+1)(c+κ+1)=0,(2.6) which can be rewritten as (μ−μ∗)2+ρ−ρ∗=0,where μ∗=1+ (c+κ)2−c 2(c+κ), ρ∗=μ∗2−κ(c+1)1+ 1 c+κ (2.7) represent the coordinates in the (μ, ρ)-plane of the vertex of the quadratic (2.6); see Figure 3. Now the application of Theorem 1.1 allows us to state the following result. Proposition 2.1. Let us consider system (2.4) and assume that c>0and the parameter κsatisfies 0<κ<κ max(c)=√c2+c−c 2.(2.8) Then the system undergoes the focus-center-limit cycle bifurcation described in Theorem 1.1 at the points of the (μ, ρ)-plane belonging to the parabolic arc defined by the quadratic equation (μ−μ∗)2+ρ−ρ∗=0(2.9) Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 3D FOCUS-CENTER-LIMIT CYCLE BIFURCATION 1939 1 κ κ MAX κ =0.25 κ 0 0.1 0.2 0.3 0.5 1 1.5 2 c Fig. 4.The graphs of the functions κmax(c)and κ1(c), which determine different regions in the plane (c, κ)as described in statements (a) and (b) of Proposition 2.1. Note the horizontal asymptote at κ=1/4. and satisfying ρ>(c+κ)μ−(c+1)κ.(2.10) The endpoints of the above parabolic arc are (μ1,ρ 1)=1−c+c 2(c+κ),c(1 −2κ)−c 2,(μ2,ρ 2)=1−c−c 2(c+κ),c(1 −2κ)+c 2, where c=c2(1 −2κ)2−4c(c+1)κ2. In the points of the above parabolic arc, the inequality D<0holds, and the following cases arise: (a) If 0<c<1and 0<κ<κ 1, where κ1=κ1(c)is the only positive root of the quartic (c+κ)4+4cκ(c+κ)−c2=0,(2.11) then μ1<μ ∗<μ 2and two subcases appear; see Figure 4. (a.1) If μ1<μ<μ ∗, then at the bifurcation points of the parabolic arc given by (2.9)–(2.10) one has γ>0. Consequently, when ρvaries, the bifurcation is supercritical and the limit cycle is orbitally asymptotically stable. (a.2) If μ∗<μ<μ 2, then γ<0at the bifurcation points of the parabolic arc. Here, when ρvaries, the bifurcation is subcritical and the limit cycle is unstable. (b) If 0<c<1and κ≥κ1,orc≥1, then all the bifurcation points of the parabolic arc (2.9)–(2.10) satisfy γ>0. Therefore, the bifurcation is supercritical and the bifurcating limit cycle is orbitally asymptotically stable. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1940 EMILIO FREIRE, ENRIQUE PONCE, AND JAVIER ROS Table 2.1 List of components for the circuit. C=C1=C2100 nF L220 mH R110 kΩ R32200 Ω R1kΩ R0220 Ω Proof. Conditions T=Tcand M>0 of Theorem 1.1 lead to MT −D= 0, which is equivalent to (2.9), and to (2.10). After some manipulations, we get the inequality (c+κ)μ2−(c+2κ)μ+(c+1)κ<0, whose discriminant, namely c2−4c2κ−4cκ2, is positive due to (2.8). In fact, this expression coincides with c2. The endpoints of the parabolic arc can be obtained by solving the equation M= 0 and (2.9). To show that D<0 at the bifurcation values, as we are working at points where MT −D= 0 along with M>0, it suffices to show that T<0, which is a trivial task. To prove statements (a) and (b), it is enough to study the sign of the coefficient γin Theorem 1. Using the condition MT −D=0,wehave γ=MT(M−m)+M(d−tM)=M[T(M−m)+d−tM], and with M>0 we get sign (γ) = sign [T(M−m)+d−tM]. Thus, using (2.5), (2.6), and canceling a factor r+μ>0, we conclude that sign (γ) = sign (c+κ)2+c+2κ−2(c+κ)μ= sign (μ∗−μ).(2.12) Assume now that 0 <c<1 and 0 <κ<κ 1. Thus, the left-hand side of (2.11) is negative, which implies (c+κ)2<c.Then μ1<μ ∗<μ 2,and statement (a) follows. When c≤1 and κ≥κ1, we have (c+κ)2≥c. If c>1, then we have (c+κ)2> c>c .In both cases, we conclude that μ∗≥μ2,and statement (b) follows. For the sake of completeness, if we define, for 0 <c<1, the constants q1=3 27c+26+3 81c2+ 156c+75>0,q 2=q1+1 q1−4>0, we obtain κ1(c)=c 6q2+c6c q2−cq2 6−2c−c, which is represented for 0 <c<1 in Figure 4. The above proposition enables us to design the electronic oscillator by choosing adequately the component values of the circuit. In particular, in order to minimize the signal distortion from the sinusoidal wave form, one must select parameters not far from the bifurcation curve where the onset of periodic oscillations has been predicted. To assess the accuracy of piecewise linear modeling for this circuit, a SPICE implementation of the circuit was made; see [QNPS93]. The values chosen for the components are in Table 2.1, while the operational amplifier used was an LM324, Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 3D FOCUS-CENTER-LIMIT CYCLE BIFURCATION 1941 M = 0 MT-D = 0 D = 0 0 0.2 0.4 0.6 0.8 1 ρ 0.2 0.4 0.6 0.8 1 μ Fig. 5.The parabolic arc (thick line) in plane (μ, ρ)predicted by Proposition 2.1 for κ=0.1. The horizontal line indicates the path followed as μvaries for the fixed value of ρused in the simulations. The dashed line represents points with D=0, so that above it we have D<0. The dotted straight line indicates points with M=0, and above it we have M>0. with a supply voltage of 9V and a measured saturation voltage of 8.5V (with slight variations around). For these values, we have c=C1 C2 =1,κ=RR0C L=0.1<√2−1 2≈0.2071, so that we can apply Proposition 2.1 and, in particular, its statement (b). Note that μ=RR2 R1R3≈R2 22000,ρ=R2C L≈0.4545, so that by varying R2we move μ, describing a horizontal path that crosses the curve corresponding to the locus of bifurcation points, as shown in Figure 5. For the above value of ρ, the bifurcation takes place for the value ¯μ≈0.4924, in accordance with (2.6), that corresponds with the value R2≈10833Ω,and oscillations will appear by increasing R2above this critical value. In Figures 6 and 7, we show the comparison between some experimental results taken from a SPICE simulation, once put into dimensionless form, and the predictions of Theorem 1.1 for the amplitude and the period of the bifurcating limit cycle. The excellent agreement achieved validates the piecewise linear model assumed for the operational amplifier nonlinear characteristic. 3. Proof of Theorem 1.1. In this section we provide the results necessary to prove Theorem 1.1. For the critical value of the bifurcation parameter Tc=D/M, the matrix AChas a pair of imaginary eigenvalues, so that for Tin a neighborhood of Tcthe eigenvalues of ACwill be α±iβ and δ∈R. The characteristic polynomial of ACis p(λ) = det (AC−λI)=−λ3+Tλ2−Mλ+D, Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1948 EMILIO FREIRE, ENRIQUE PONCE, AND JAVIER ROS Due to the symmetry of the orbit, its period is equal to 2(τC+τL). Substituting expansion (3.18) into (3.9), and computing the above expression for the period, we get the expansion given for Per. We will now determine the amplitude of the periodic orbit. By using the variation of parameters formula, the solution of system (1.2) in zone Ris x(τ)=eALτx3(τL)+τ 0 eAL(τ−s)b(τL)ds,(3.19) so that its first component is x1(τ)=eT 1⎧ ⎨ ⎩eALτ⎡ ⎣ 1 −x1 2(τL) −x1 3(τL)⎤ ⎦+ ∞ ! i=0 Ai L τi+1 (i+ 1)!"b(τL)⎫ ⎬ ⎭.(3.20) Let τ∗be the time when |x1|attains its maximum value in zone R. Taking derivatives with respect to τin (3.20), and imposing that it must vanish at τ∗, we get G(τL,τ∗)= dx1(τ) dτ τ=τ∗ =eT 1eALτ∗⎡ ⎣ x1 2(τL)+T(τL) x1 3(τL)+M D⎤ ⎦=0.(3.21) Now using expressions (3.11) and (3.13) and computing the power series of Gin (τL,τ∗)at(0,0), we obtain G(τL,τ∗)=M 2τL−Mτ∗+D−Mt 12 τ2 L+Mt−D 2τLτ∗+D−Mt 2τ∗2+O(τL,τ∗)3. Hence, (3.21) defines implicitly in a neighborhood of (0,0) a function τ∗=ψ(τL) such that G(τL,ψ(τL)) = 0, namely, τ∗=1 2τL+Mt−D 24Mτ2 L+O(τ4 L). Substituting the above expansion together with (3.11), (3.13), and (3.17) into the expression (3.20), we get a=x1(τ∗)=1+M 8τ2 L+1 1152M(13D2−11DMt +15M2m−2M2t2)τ4 L+O(τ5 L). Using expression (3.18) for τL, we obtain the final expression for the amplitude a. Let us now compute the characteristic multipliers of the bifurcating limit cycle. Due to the similarity relationship established in Proposition 3.2, we conclude that the product exp(ALτL)·exp(ACτC) corresponding to a solution of (3.5) has an eigenvalue equal to −1. We will denote by λrand λathe other two eigenvalues that correspond to the eigenvalues of the derivative DpπLC (p0) of the transition map associated with the semiorbit. The product of the three eigenvalues is then equal to −λrλa= det eALτLdet eACτC. Using that det(eAτ ) = exp(τtrace(A)),we get −λrλa=eτLt+τCT.(3.22) Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 3D FOCUS-CENTER-LIMIT CYCLE BIFURCATION 1949 The expansion of the product of exponentials in (3.16) leads to an expression of the form eALτLeAC(τL)τC(τL)=H0+τLH1+τ2 LH2+···.(3.23) To compute the above matrices Hi, we write I+ALτL+A2 L τ2 L 2! +··· × eAC(0)τC(0) +τL d dτL eAC(τL)τC(τL)τL=0 +τ2 L 2! d2 dτ2 L eAC(τL)τC(τL)τL=0 +···". From expansions (3.6)–(3.13), we obtain τC(0) = π/M1/2,τ C(0) = −1, τ C(0)=0, A C(0) = A C(0) = 0,and using these values in the above expression, we finally get H0=eAC(0)τC(0) =⎡ ⎣ D2K−1−DMK M2K 0−10 D2MK −DM2KM 3K−1⎤ ⎦, H1=⎡ ⎣ t−D/M m−M d−D⎤ ⎦D2K−1−DMK M2K and H2=Mt−D 2MH1, where Khas been defined in (3.15), and it is emphasized that H1and H2are rank-one matrices. The matrix H0has eigenvalues −1 (double) and λ0= exp(πD/M3/2). For the single eigenvalue λ0, we select a right eigenvector v0=[1,0,M]Tand a left eigenvector wT 0=[D2/M 2,−D/M, 1]. We will denote by λathe eigenvalue of DπLC (p0) that for τL= 0 is equal to λ0. Since the eigenvalue λ0of H0is simple, we can apply perturbation theory (see section 2.8 of [Wi65]) to assure that the equality H0+τLH1+τ2 LH2+···v0+τLv1+τ2 Lv2+··· =(λ0+τLλ1+τ2 Lλ2+···)v0+τLv1+τ2 Lv2+··· holds for certain vectors v1,v2... . As a consequence of Proposition 3.2 and (3.23), we get λa=λ0+τLλ1+τ2 Lλ2+···. After some computations, we arrive at λ1=wT 0H1v0 wT 0v0 =Mt−D M+γM D2+M3λ0, λ2=wT 0(H2v0+(H1−λ1I)v1) wT 0v0 =(Mt−D)D2+M3+γ 2M2(D2+M3)2(λ0+1) ×(Mt−D)D2+M3λ0+D2−M3+2M2(dM −Dm)λ0. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1950 EMILIO FREIRE, ENRIQUE PONCE, AND JAVIER ROS The logarithms μrand μaof characteristic multipliers of the complete periodic orbit must satisfy eμr=λ2 r,e μa=λ2 a,(3.24) while from (3.22) we get the relationship μr+μa=2tτL+2TτC.(3.25) From (3.24) and using the computed simple eigenvalue λa, we obtain μa= 2 log λ0+λ1τL+λ2τ2 L+O(τ3 L)= 2 log λ0+ 2 log 1+λ1 λ0 τL+λ2 λ0 τ2 L+O(τ3 L) =2λ0+2λ1 λ0 τL+2λ2 λ0−λ2 1 λ2 0τ2 L+O(τ3 L). Substituting here λ1and λ2, and using expansion (3.18) of τL, we finally get the expression for μathat appears in Theorem 1.1. Using in (3.25) the expansions (3.9) for τCand (3.18) for τL, we compute μr. Since the last assertion of Theorem 1.1 is a direct consequence of previous statements, its proof is now completed. Acknowledgments. The authors sincerely appreciate the careful reading of the anonymous referees and their interesting suggestions that have notably improved the paper. They also want to acknowledge Jorge Gal´an for his invaluable suggestions and comments on a preliminary version of the manuscript, Manuel Rom´an for his help with SPICE simulation of the circuit, and Fernando Fern´andez for helping with figures. REFERENCES [AVK66] A. A. Andronov, A. A. Vitt, and S. E. Khaikin,Theory of Oscillators, Dover, New York, 1966. [CFPT02] V. Carmona, E. Freire, E. Ponce, and F. Torres,On simplifying and classifying piecewise-linear systems, IEEE Trans. Circuits Systems I Fund. Theory Appl., 49 (2002), pp. 609–620. [CH82] S. N. Chow and J. K. Hale,Methods of Bifurcation Theory, Springer-Verlag, New York, Berlin, 1982. [CWHZ93] L. Chua, C. Wu, A. Huang, and G. Zhong,A universal circuit for studying and generating chaos—Part I: Routes to chaos, IEEE Trans. Circuits Systems I Fund. Theory Appl., 40 (1993), pp. 732–744. [FGA84] E. Freire, L. Garc´ ıa-Franquelo, and J. Aracil,Periodicity and chaos in an autonomous electronic oscillator, IEEE Trans. Circuits Systems, 31 (1984), pp. 237–247. [FRGP93] E. Freire, A. J. Rodr´ ıguez-Luis, E. Gamero, and E. Ponce,A case of study for homoclinic chaos in an autonomous electronic circuit: A trip from Takens– Bogdanov to Hopf–˘ Sil’nikov, Phys. D, 62 (1993), pp. 230–253. [FPR99] E. Freire, E. Ponce, and J. Ros,Limit cycle bifurcation from center in symmetric piecewise-linear systems, Internat. J. Bifur. Chaos, 9 (1999), pp. 895– 907. [GK92] M. G. M. Gomes and G. P. King,Bistable chaos. II. Bifurcation analysis, Phys. Rev. A, 46 (1992), pp. 3100–3110. [HBCJM91] J. J. Healey, D. S. Broomhead, K. A. Cliffe, R. Jones, and T. Mullin,The origins of chaos in a modified Van der Pol oscillator, Phys. D, 48 (1991), pp. 322–339. [Ke93] M. P. Kennedy,Three steps to chaos—Part II: A Chua’s circuit primer, IEEE Trans. Circuits Systems I Fund. Theory Appl., 40 (1993), pp. 657–674. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 3D FOCUS-CENTER-LIMIT CYCLE BIFURCATION 1951 [Kr87] G. A. Kriegsmann,The rapid bifurcation of the Wien bridge oscillator, IEEE Trans. Circuits Systems, 34 (1987), pp. 1093–1096. [Ma93] R. Madan,Chua’s Circuit: Paradigm for Chaos, World Scientific, Singapore, 1993. [MGHLVM03] M. B. Monagan, K. O. Geddes, K. M. Heal, G. Labahn, S. M. Vorkoetter, J. McCarron, and P. DeMarco,Maple 9Introductory Programming Guide, Maplesoft, Waterloo, ON, 2003. [QNPS93] T. Quarles, A. R. Newton, D. O. Pederson, and A. SangiovanniVincentelli,Spice3Version 3f3 User’s Manual, Department of Electrical Engineering and Computer Sciences, University of California Berkeley, Berkeley, CA, 1993. [Ro03] J. Ros,Estudio del Comportamiento Din´amico de Sistemas Aut´onomos Tridimensionales Lineales a Trozos, Ph.D. dissertation, Universidad de Sevilla, Seville, Spain, 2003 (in Spanish). [SYM81] R. Shinriki, M. Yamamoto, and S. Mori,Multimode oscillations in a modified Van der Pol oscillator containing a positive nonlinear conductance, IEEE Proc., 69 (1981), pp. 394–395. [Wi65] J. H. Wilkinson,The Algebraic Eigenvalue Problem, Oxford University Press, Oxford, UK, 1965. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php