Bifurcation set for a disregarded Bogdanov-Takens unfolding: Application to 3D cubic memristor oscillators

arXiv:1811.04862v1 [math.DS] 12 Nov 2018 1 BIFURCATION SET FOR A DISREGARDED BOGDANOV-TAKENS UNFOLDING. APPLICATION TO 3D CUBIC MEMRISTOR OSCILLATORS. ANDR´ ES AMADOR1AND EMILIO FREIRE2AND ENRIQUE PONCE3 Abstract. We derive the bifurcation set for a not previously considered threeparametric Bogdanov-Takens unfolding, showing that it is possible express its vector field as two different perturbed cubic Hamiltonians. By using several firstorder Melnikov functions, we obtain for the first time analytical approximations for the bifurcation curves corresponding to homoclinic and heteroclinic connections, which along with the curves associated to local bifurcations organize the parametric regions with different structures of periodic orbits. As an application of these results, we study a family of 3D memristor oscillators, for which the characteristic function of the memristor is a cubic polynomial. We show that these systems have an infinity number of invariant manifolds, and by adding one parameter that stratifies the 3D dynamics of the family, it is shown that the dynamics in each stratum is topologically equivalent to a representant of the above unfolding. Also, based upon the bifurcation set obtained, we show the existence of closed surfaces in the 3D state space which are foliated by periodic orbits. Finally, we clarify some misconceptions that arise from the numerical simulations of these systems, emphasizing the important role played by the existence of invariant manifolds. Bifurcation set, Bogdanov-Takens, homoclinic orbit, heteroclinic connection, Melnikov function, Memristor oscillators 1. Introduction In planar systems, the existence of some local bifurcations may reveal the presence of other bifurcations of global character [13, 14] and the curves that determine these global phenomena are difficult to determine. This is, for instance, the case regarding the appearance of homoclinic or heteroclinic connections. A homoclinic connection is an orbit of the system that joins a saddle equilibrium point to itself, and generally creates or destroys periodic orbits (see, for instance [35]). A heteroclinic connection joins two different equilibrium points of a system and the existence of this connection can determine changes in the basin of attraction of a positively invariant set. Following [18], the techniques to study homoclinic orbits in planar vector fields were well developed during the 1920s in the works of Dulac. The fundamental idea is that the recurrent behavior near a connecting orbit should be studied in a fashion similar to that used in studying periodic orbits via a Poincar´e return map. But there are some additional complications in the study of homoclinic orbits compared to that of periodic orbits which significantly complicate the analysis. 1Facultad de Ingenier´ıa y Ciencias, Departamento de Ciencias Naturales y Matem´aticas, Pontificia Universidad Javeriana-Cali, Cali, Colombia. 2,3Departamento de Matem´atica Aplicada, Escuela T´ecnica Superior de Ingenier´ıa, Avda. de los Descubrimientos, 41092 Sevilla, Spain. [email protected], [email protected] 3eponce[email protected] 1 2 A. AMADOR AND E. FREIRE AND E. PONCE Usually, the bifurcation curves of homoclinic and heteroclinic connections are studied by numerical continuation techniques [17, 18, 9, 36]. On the other hand, when a planar system can be written as a perturbed Hamiltonian system, we can calculate some Melnikov functions, introduced by Melnikov in [28], and under certain hypotheses, the zeros of the associated Melnikov function determine the existence of periodic orbits, homoclinic loops or heteroclinic connections, see for instance [8, 20, 3]. As show later, we will resort to such Melnikov functions for setting information on such global bifurcations curves in a two-parametric plane for a specific family of differential systems. The normal form of the Bogdanov-Takens (see [35]) bifurcation is given by ˙x=y, ˙y=µ1+µ2x+x2±xy. Following the classification proposed in [14], the deformation of codimension three of the previous normal form is given by the unfolding ˙x=y, ˙y=µ1+µ2x+αx3+y(µ3+µ4x±x2), (1) The presence of this type of systems has been reported in different applications, see [16, 2, 24]. On the study of bifurcation phenomena in these systems many contributions have been made. In [15], the authors studied the global bifurcation diagram of the three-parameter family ˙x=y, ˙y=µ1+µ2x−x3+y(µ3−3x2), and fixing µ3>0, they obtained analytical approximations to the bifurcation curves of the homoclinic orbits, by using Melnikov functions. Later on, the above work was quoted in [23], where a numerical analysis of the same model was performed. In [10, 12, 11], it is considered the system ˙x=y, ˙y=µ1+µ2x−x3+y(µ3+µ4x−x2), (2) and the authors showed that it can be written as a perturbed Hamiltonian system, reporting the maximum number of limit cycles. Later, by taking the parameter µ4= 0 in (2), the authors in [4, 5, 7, 6] analyzed the system as a Li´enard system, its local bifurcations were characterized, and a numerical study of the global bifurcations was done. While all the above references dealt with the focus case, in this work we study the saddle case ˙x=y, ˙y=µ1+µ2x+x3+y(µ3−3x2), (3) which up to the best of our knowledge, seems to be a disregarded case with rather interesting dynamic behavior. In fact, our motivation comes from the analysis of certain 3D memristor oscillators [1, 30, 22], where under specific hypotheses on the memristor characteristics, such system appears in a natural way after a dimensional reduction achieved thanks to the existence of a first integral. BIFURCATION SET FOR A DISREGARDED BOGDANOV-TAKENS UNFOLDING 3 The paper is organized in the following way. First, in section 2 we review the information that can be gained by means of a local analysis of the system. Our main results appear in section 3, where we apply Melnikov theory to approximate the homoclinic and heteroclinic curves on a convenient parameter plane. We obtain the global bifurcation set for system (3), by splitting such a plane in regions with different qualitative dynamical behavior. Next, in section 4, We show how the above analysis is useful for deriving all the possible responses of certain 3D canonical memristor oscillators, when the flux-charge characteristics function is a specific cubic polynomial, generalizing some results given in [1]. As one of the possible dynamical behaviors, we focus or attention in showing, thanks to the previous analysis, the existence of a topological sphere in the 3D phase-space completely foliated by periodic orbits. Thus, we confirm previous numerical results reported in [29, 25]. The necessity of incorporating rigorous techniques in the analysis of memristor oscillators is emphasized with the material of section 5, where following a similar procedure for the dimensional reduction of section 4, we can refute several recently published studies that report the existence of an infinite number of hidden attractors in a three-dimensional memristor-based autonomous Duffing oscillator. Some technical results are relegated to the appendix. 2. Local bifurcations In this section, we study the local bifurcations that occur in system (3). First, we note that the system is invariant under the transformation (x, y, µ1, µ2, µ3)→(−x, −y, −µ1, µ2, µ3). Therefore, it is sufficient to study the bifurcation diagram for µ1>0.The equilibrium points of the system are of the form (x, y) = (˜x, 0), being ˜xa solution of the cubic µ1+µ2x+x3= 0,and its Jacobian matrix of is given by (4) J(x, y) = 0 1 µ2+ 3x2−6yx µ3−3x2. Remark 1. Note that for µ3≤0the divergence of system (3) does not change sign, thus from Bendixson’s criterion [21], the system does not have periodic solutions. First, we provide a technical result that provides a study of the number of equilibria in system (3) and their topological nature. Lemma 2. Consider system (3), the following statements hold. (a) If µ2≥0or we have µ2<0with 27µ2 1+ 4µ3 2>0,then the system has only one equilibrium point. (b) If µ2<0and 27µ2 1+ 4µ3 2= 0 we have two equilibrium points. (c) If µ2<0and 27µ2 1+ 4µ3 2<0,then the system has three equilibrium points xi= (si,0) with i∈ {L, C, R}such that (5) sL<−(−µ2/3)1/2< sC<(−µ2/3)1/2< sR, and sL+sC+sR= 0.Furthermore, xLand xRare saddles while xCis an antisaddle (node or focus). Proof. We study the roots of the polynomial p(x) = µ1+µ2x+x3.Since p′(x) = µ2+ 3x2,if µ2≥0 we obtain p′(x)≥0 and so the polynomial has only one root. 4 A. AMADOR AND E. FREIRE AND E. PONCE In the rest of the proof we assume µ2<0. The derivative p′(x) vanishes at the points x±=±(−µ2/3)1/2being a maximum and minimum local respectively, also a direct computation gives p(x±) = µ1∓2(−µ2/3)3/2. When p(x−)<0 or p(x+)>0 the graph of p(x) only crosses once the x-axis and so these inequalities provides the condition 27µ2 1+ 4µ3 2>0,and the statement (a) follows. Assuming p(x−) = 0 or p(x+) = 0,the statement (b) follows. Finally, if p(x+)<0< p(x−) then we have three roots as indicated in (5). Moreover, using the relation between roots and coefficients of polynomials, we get µ1=−sLsCsR, µ2=sCsL+sCsR+sLsR, sL+sC+sR= 0. For the Jacobian matrix given in (4), we get J(si,0) = −(µ2+ 3s2 i) = −p′(si).As we know that p′(sL)>0, p′(sC)<0 and p′(sR)>0,the conclusion follows and the proof is complete.  In the next result, we give a characterization of the local bifurcations of system (3) on the parametric plane (µ2, µ1),assuming a fixed value of the parameter µ3. Note that wee have put the µ1axis in the plane (µ2, µ1) vertically, being the µ2axis the horizontal one. Proposition 3. The following statements hold for system (3). (a) Given µ3∈Rthe parameter values in the set (6) ϕsn ={(µ2, µ1) : 27µ2 1+ 4µ3 2= 0}, correspond with saddle-node bifurcation points of equilibria. In particular, the system has a cusp bifurcation of equilibria at µ2=µ1= 0. (b) Given µ3>0,the parameter values in the set (7) ϕH={(µ2, µ1) : µ1=±(µ3/3)3/2∓(µ3/3)1/2µ2, µ2<−µ3}, represent Andronov-Hopf bifurcation points of codimension one for the central equilibrium point xC,see Lemma 2(c). (c) The set defined in (7) determines a symmetric pair of straight half lines emanating from two points corresponding to Bogdanov-Takens bifurcation points, namely (8) BT±≡−µ3,±2 (µ3/3)3/2. Proof. Statement (a) is a direct consequence of the equations x3+µ2x+µ1= 0, µ2+ 3x2= 0, to be fulfilled for any non-hyperbolic equilibrium (x, 0) at a saddle-node bifurcation. Let (˜x, 0) an equilibrium point of system (3). Considering the Jacobian matrix J given in (4), then J(˜x, 0) has two purely imaginary eigenvalues when taking µ3>0, the value ˜xsatisfies ˜x=±pµ3/3 with µ2<−µ3<0,because then µ2+ 3˜x2<0. The last inequality is fulfilled only for the equilibrium point xC,see Lemma 2(c). Since (˜x, 0) is an equilibrium point we have µ1+µ2±pµ3/3+±pµ3/33= 0, and statement (b) follows. To show statement (c) is sufficient to consider the equations trace (J(˜x, 0)) = det(J(˜x, 0)) = 0. BIFURCATION SET FOR A DISREGARDED BOGDANOV-TAKENS UNFOLDING 5 3. Global bifurcations In this section we will complete the bifurcation analysis of system (3). We will write the system as a perturbed Hamiltonian to which the Melnikov theory can be applied. This can be done in different ways, as indicated in the next result. The possibility of resorting to one of the two next reparametrization forms will be helpful later. Proposition 4. System (3) can be written as two different perturbed Hamiltonian systems, as follows. (a) Taking (9) µ1=ε4ν1, µ2=−ε2ν2, µ3=ε2ν3, the system can be rewritten as ˙x=y, ˙y=−ν2x+x3+εν1+ν3y−3x2y, (10) which for ε= 0 corresponds to the Hamiltonian (11) H1(x, y) = y2 2+ν2 x2 2−x4 4. (b) Taking (12) µ1=ε3ν1, µ2=−ε2ν2, µ3=ε2ν3, the system can be rewritten as ˙x=y, ˙y=ν1−ν2x+x3+ε(ν3y−3x2y), (13) which for ε= 0 corresponds to the Hamiltonian (14) H2(x, y) = y2 2−ν1x+ν2 x2 2−x4 4. Proof. The blow-up transformation x1= (1/ε)x, y1= (1/ε2)y, and ˜ t=εt, allows to rewrite system (3) as x′ 1=y1, y′ 1=x3 1+µ2 ε2x1+µ1 ε3+µ3 εy1−3εx2 1y1, where the prime denotes derivatives with respect to the new time ˜ t. Now, using (9) and (12), after some elementary algebra we obtain systems (10) and (13), respectively.  The phase portrait for the unperturbed Hamiltonian systems (10) and (14) are shown in Figure 1. Note that when ν1= 0 we obtain H1(x, y) = H2(x, y),and so in that case it is sufficient to study the properties of the Hamiltonian H1. Now, we will consider the heteroclinic connections of unperturbed Hamiltonian system (11). The Hamiltonian has a pair of heteroclinic connections Γ±(t) = (x(t),±y(t)),parameterized by x(t) = √ν2tanh pν2/2t, y(t) = ν2 √2sech2pν2/2t, (15) 6 A. AMADOR AND E. FREIRE AND E. PONCE where −∞ < t < ∞and ν2>0.In the next result, we compute the Melnikov function along the heteroclinic connection Γ+for the unperturbed Hamiltonian system (10), and by using (9), we obtain the approximate bifurcation curves for heteroclinic connections for system (3). Figure 1. (a) Phase portrait of unperturbed Hamiltonian system (10) with ν2= 0.2. We show in green the two heteroclinic orbits, while the non-closing stable and unstable manifolds of the saddle points are shown in red. (b) Phase portrait of unperturbed Hamiltonian system (13) with ν1= 0.3 and ν2= 1. We draw in green the homoclinic orbit, the stable and unstable manifolds of the saddle points are shown in red. Proposition 5. If we consider perturbed Hamiltonian system (10) and ν= (ν1, ν2, ν3) with ν2>0and ν2 2/4< ν3 2/27 (see Lemma 2(c)) then the Melnikov function along of the heteroclinic connection Γ±is given by (16) Mht(ν) = 2 15√ν215ν1+ 5√2ν2ν3−3√2ν2 2. Proof. The system can be written as ( ˙x, ˙y)T=f(x, y) + εg(x, y), where f(x, y) = (y, −ν2x+x3)Tand g(x, y) = (0, ν1+ν3y−3x2y)T.Thus, we have f∧g=y(ν1+ν3y−3x2y). Accordingly, the Melnikov function is defined by Mht(ν) = Z∞ −∞ f(x(t),±y(t)) ∧g(x(t),±y(t))dt = =Z∞ −∞ ±y(t)ν1±(ν3−3x2(t))y(t)dt, where x(t) and y(t) are defined as in (15). After a direct computation we obtain (16).  BIFURCATION SET FOR A DISREGARDED BOGDANOV-TAKENS UNFOLDING 7 By using the Melnikov theory and by fixing one parameter of system (3), we can give an approximation of the heteroclinic connection curves in the remaining parameters plane. Proposition 6. Consider system (3) with µ3>0sufficiently small and the parametric plane (µ2, µ1). Then the system has a unique hyperbolic heteroclinic connection in a neighborhood of the curve (17) ϕht ={(µ2, µ1)∈R2:µ1=±√2 15 µ2(3µ2+ 5µ3), µ26=−5µ3/3}. Proof. Fixing ν3= 1 in the Melnikov function given in (16), and imposing the condition Mht(ν1, ν2) = 0,we obtain ν1=√2 15 ν2(3ν2−5) . From (9) we get ε=√µ3,µ1=µ2 3ν1,and µ2=−µ3ν2,so that µ1=µ2 3ν1=−µ2 3 √2 15 µ2 µ33−µ2 µ3−5, and the conclusion follows.  When µ1= 0 and µ3>0 on the parameter plane (µ2, µ1),we obtain the point of double heteroclinic connections (18) DHT ≡(−5µ3/3,0) . We recall that Schecter’s points are co-dimension two points defined by the intersection of a saddle-node curve and a homoclinic or heteroclinic curve, for more details see [32]. Taking the intersection points of the saddle-node bifurcation curve and the heteroclinic curves given in (6) and (17) respectively, we obtain a first-order approximation of Schecter’s points of the system. Since the system is symmetric with respect to the parameter µ1,the system has four Schecter’s points (see Figure 3), these points are S± 1≡ρ1(5/27),∓(5√10/729) p18µ3+ 5 + √5, S± 2≡ρ2(5/27),±(5√10/729) p18µ3+ 5 −√5, (19) where ρ1=9µ3+ 5 −√5p18µ3+ 5, ρ2=√5p18µ3+ 5 −9µ3−5. Now, by using the homoclinic connection of Hamiltonian system (14), we compute the associated Melnikov function for system (3) when ν3= 1. Proposition 7. If we consider system (13) and ν= (ν1, ν2, ν3)with ν1>0,ν2>0 and ν3= 1, then the Melnikov function associated to the homoclinic orbit with connection point (0, sR), it is given by (20) M(ν) = √2cosh2(θ) cosh2(θ) + 2 (F1(θ) + ν2F2(θ)) , 8 A. AMADOR AND E. FREIRE AND E. PONCE where F1(θ) =720θ−320 sinh θ+ 240θcosh3θ−320 cosh2θsinh θ− −80 cosh4θsinh θ+ 480θcosh θ, F2(θ) =1440θcosh θ−768 sinh θ−cosh3θ−1344 cosh2θsinh θ− −48 cosh4θsinh θ (21) and 0< θ < ∞, with cosh θ=2s ω, ω2= 2(ν2−s2 R)>0, ν1=ν2sR−s3 R, being sRthe biggest positive root of the equation ν1−ν2x+x3= 0, see Figure 2. Proof. We consider the unperturbed Hamiltonian system given in (13) with ν2>0. From Lemma 2(c), the system has 3 equilibrium points xi= (si,0), where xLand xRare saddle points and xCis a focus or node and sL< sC< sR, sL+sC+sR= 0, sLsCsR=−ν1. We study only the case ν1>0, for the case ν1<0 is analogous. System (13) can written as ( ˙x, ˙y)T=f(x, y) + εg(x, y). Now, assuming ν1>0,by Green’s Theorem, the homoclinic Melnikov function of the system can rewritten as Mh(ν) = Z ZD(ν1,ν2)−∂g(x, y) ∂y dA, where Dis the region bounded by the homoclinic orbit which joins the equilibrium point (sR,0) to itself. By fixing ν3= 1 (that is µ3>0), and taking p(x) = ν1−ν2x+x3,we get p(sR) = ν1−ν2sR+s3 R= 0,that is (22) ν1=sR(ν2−s2 R), and so ν2−s2 R>0.Taking the auxiliary function q(x) = Zx 0 p(x)dx =ν1x−ν2 x2 2+x4 4, and using (14), the homoclinic loop is given by the points (x, y± s(x)) where x≤x≤ sR, y± s(x) = ±√2pq(x)−q(sR), and y± s(x) = y± s(sR) = 0,see Figure 2. Now, the Melnikov function is thanks to the symmetry of the loop Mh(ν) = 2 ZsR x (3x2−1)dx Zy+ s(x) 0 dy = =√2Zs x (3x2−1)(sR−x)q(x+sR)2−2(ν2−s2 R)dx = =√2Zs x (3x2−1)(sR−x)p(x+sR)2−ω2dx, BIFURCATION SET FOR A DISREGARDED BOGDANOV-TAKENS UNFOLDING 9 where from (22) ω2= 2(ν2−s2 R),and we have used that q(x)−q(sR) = 1 4(x−sR)2(x+sR)2−ω2. Taking the change of variable x+sR=ωcosh θ, and noting that q(sR)−q(x) = 0 we see that x+sR=ω, which corresponds to θ= 0, while for x=sRthe corresponding values of θ=θsRsatisfy cosh θsR= 2sR/ω, or also ω2cosh2θsR= 4s2 R,that is (s2 R+ν2) cosh2θsR= 2s2 R,and so we get s2 R=cosh2θsR 2 + cosh2θsR ν2. Now we arrived to Mh(ν) = √2ω2ZθsR 0 (1 −3(ωcosh θ−sR)2)(2sR−ωcosh θ) sinh2θdθ, and after some computations, we obtain (20) and (21), where θsRhas been simplified to θ. 0 x 0 y y+ s(x) y− s(x) sC sLxsR Figure 2. Homoclinic orbit which joins the saddle equilibrium point (sR,0) to itself. As a direct consequence of the above result, we give an analytical approximation of the bifurcation curves for homoclinic connections of system (3). Proposition 8. Consider system (3) with µ3>0sufficiently small and the parametric plane (µ2, µ1).Then the system has a unique homoclinic orbit in a neighborhood of the curve (23) ϕh={(µ2, µ1)∈R2:µ2=−µ3ν2(θ), µ1=±µ3/2 3ν1(θ),0< θ < ∞}, where ν2(θ) = 10(cosh 2θ+ 5)(9 sinh θ+ sinh 3θ−12θcosh θ) 3(370 sinh θ+ 115 sinh 3θ+ sinh 5θ−60θ(11 cosh θ+ cosh 3θ)), ν1(θ) = ν2(θ)s−s3, s2=cosh2θ 2 + cosh2θν2(θ). (24) 16 A. AMADOR AND E. FREIRE AND E. PONCE system (46). This result guarantees the existence of a topological sphere in the 3D phase-space completely foliated by periodic orbits. Proposition 11. Consider system (30) with β, ξ > 0,, the function qdefined as in (31),a2−3b < 0and (37) a2−3b+ 3β > 0 sufficiently small. Additionally, suppose that the following inequalities hold 0<3b−a2<3ξ/β, β(3b−a2)−3ξ < (5/2) 3b−a2−3ββ < 0. (38) Then for all h∈Rwith −A 27 < h < B 27, where A=4a2β+ 3β2−12bβ + 9ξpa2−3b+ 3β+ 9aξ + 2a3β−9abβ, B=4a2β+ 3β2−12bβ + 9ξpa2−3b+ 3β−9aξ −2a3β+ 9abβ, (39) the system has a stable periodic orbit. Moreover, there exist a topological sphere Ω (see Figure 10) foliated by such periodic orbits. Proof. From Remark 10, and after substituting the values of µ1, µ2and µ3given in (35), we obtain the inequalities (37)-(38). Now, from (36) we obtain |µ1|< (1/3) (µ3/3)1/2(µ3−3µ2),so that µ3−3µ2>0,since from hypotheses we have µ2<−(5/2)µ3<0.Now after some algebra we obtain |27h+ 9aξ + 2a3β−9abβ|<4a2β+ 3β2−12bβ + 9ξpa2−3b+ 3β. Taking into account the absolute value, and grouping terms, we obtain the values of Aand Bdefined in (39). Finally, from Remark 10 system (30) has a stable periodic orbit on each Shdefined in (32), so varying the parameter h, we obtain a sphere foliated by such periodic orbits.  5. False Hidden Attractors in Memristor-Based Autonomous Duffing Oscillators An attractor is called a hidden attractor if its basin of attraction does not intersect any neighborhood of equilibria; otherwise, it is called a self-excited attractor, for more details see [26, 27]. Recently in [33, 34, 19] it was reported the existence of an infinite number of hidden attractors in a memristor-based autonomous Duffing oscillators, whose memristance function is a cubic polynomial. Here, by using a similar approach to the followed in the previous section, we will show that such hidden attractors are not possible, so that the numerical simulations included in [33, 34, 19] are misleading. The quoted memristor based autonomous Duffing oscillator is modeled by the dynamical system ˙x=y, ˙y=z, ˙z=−αz −M(x)y, (40) BIFURCATION SET FOR A DISREGARDED BOGDANOV-TAKENS UNFOLDING 17 Figure 10. Using (62) on each invariant manifold Shdefined in (32), some slices of the surface Ω given by Proposition 11 for system (30) with parameters a= 1, b = 4.8, β = 5 and ξ= 80 are shown. For this set of parameters we get µ3= 0.106 >0, µ2=−2.3<−(5/2)µ3, A= 1180.1 and B= 152.2. where the memristance function M(possibly discontinuous) is defined as (41) M(x) = dφ(x) dx and φis a continuous function. System (40) has a continuum of equilibria, since any point of the x-axis is an equilibrium point. In the next result, we show that even system (40) does not belong to the family (46) of the appendix, the system also has the property of possessing an infinite number of invariant manifolds. Proposition 12. Consider system (40) with the function Mdefined as in (41). For any h∈Rthe set (42) Sh={(x, y, z)∈R3:H(x, y, z) = h} is an invariant manifold for the system, where we have introduced the continuous function (43) H(x, y, z) = φ(x) + αy +z. Therefore, the system has an infinite number of invariant manifolds foliating the whole R3, and so the dynamics is essentially two-dimensional. Proof. Taking Has in (43), define for any solution (x(τ), y(τ), z(τ)) of (40) the auxiliary continuous function h(τ) = H(x(τ), y(τ), z(τ)) Now, a direct computation gives, excepting the points of possible non-differentiability, h′(τ) = dφ(x) dx ˙x+α˙y+ ˙z=M(x)y+αz −αz −M(x)y= 0. Then his piecewise constant along the orbits of (40), but as his continuous by definition, it should be globally constant. In short, the level sets of Hare invariant for the flow.  18 A. AMADOR AND E. FREIRE AND E. PONCE Now, by using the above result, we reduce the study of the dynamical behavior of the system, to the study of a planar system. Proposition 13. Consider system (40) with the function Mdefined as in (41). Then on each invariant set Shdefined in (42) the system is topologically equivalent to the planar system ˙x=y, ˙y=−φ(x)−αy +h. (44) Moreover, (x(τ), y (τ)) ∈R2is a solution of the above system if and only if Eh(x(τ), y (τ)) is a solution of system (40), where (45) Eh(X(τ), Y (τ)) =  x(τ) y(τ) h−φ(x(τ)) −αy(τ)   Proof. From Proposition 12 we can solve for zin the equation H(x, y, z) = h, and write z=h−φ(x)−αy. Replacing this expression into the first and second equation of (40) we obtain system (44). Suppose that (x(τ), y (τ)) ∈R2is a solution of system (44). Taking z(τ) = h−αy(τ)−φ(x(τ)) we obtain ˙z(τ) = −α˙y(τ)−dφ(x(τ)) dx ˙x(τ) = −α(h−φ(x(τ)) −αy(τ)) −M(x(τ))y(τ) = =−α(z(τ)) −M(x(τ))y(τ). and the proposition follows.  In the following result, we show that for α6= 0, the system does not have periodic solutions. Proposition 14. Consider system (44). The following statements hold. (a) For α= 0 the system is Hamiltonian. (b) For α6= 0 the system does not have periodic solutions. Proof. The divergence of the system is ∆ = −α. Then, when α= 0 the system corresponds to the Hamiltonian H(x, y) = y2 2+φ′(x). For α6= 0 the divergence of system (44) does not change sign, thus from Bendixson’s criterion [21], system (40) does not have periodic solutions.  Remark 15. Note that as a consequence of proposition 13 and 14, the 3Dsystem (40) system cannot have periodic orbits for any continuous function φand α6= 0. However, when α= 0 the system could have an infinite number of periodic orbits on each invariant set Shdefined in (42). BIFURCATION SET FOR A DISREGARDED BOGDANOV-TAKENS UNFOLDING 19 Regarding [33, 34, 19] authors consider system (40) with the function φ(x) = ωx +βx3and the set of parameters α= 0.0001, ω = 0.35, β = 0.85.In both quoted references, authors reported the existence of an infinite number of stable periodic orbits coexisting with an infinite number of stable equilibria, by taking into account several numerical simulations, so concluding the existence of hidden attractors. From Propositions 12 and 13, we obtain the invariant manifolds Sh={(x, y, z)∈R3:ωx +βx3+αy +z=h}, and the planar system ruling the dynamics on each Shgiven by ˙x=y, ˙y=−ωx −βx3−αy +h. From Remark 15 we note that, the system cannot have periodic orbits, and so, the statement made in the quoted papers is clearly wrong, probably after giving too much credit to numerical simulations. This emphasizes the relevance of the approach followed in this work which allows to avoid misconceptions coming just from numerical simulations. 6. Conclusions Motivated by the dynamical analysis of 3D memristor oscillators whose nonlinear characteristics is a cubic polynomial, and after showing that their dynamics is essentially two-dimensional, the need to consider a disregarded unfolding of the Bogdanov-Takens singularity naturally arose. The corresponding bifurcation set, including both local and global bifurcations has been described. While local bifurcations can be easily detected, the characterization of global bifurcations parameters curves is much more involved; only by resorting to Melnikov’s theory it was possible to obtain such curves providing a complete description of the bifurcation set. Regarding the considered 3D memristor oscillators, and by working within some parameters regions of the above bifurcation set, it has been possible to show rigorously the existence of multiple periodic orbits leading to a topological sphere. When the same approach is applied to a different family of 3D memristor oscillators, it has been shown that the oscillations are not possible, contrarily to what had been recently claimed. Acknowledgements The first author is supported by Pontificia Universidad Javeriana Cali-Colombia. E. Freire and E. Ponce are partially supported by MINECO/FEDER grant MTM2015- 65608-P and by the Consejer´ıa de Econom´ıa, Innovaci´on, Ciencia y Empleo de la Junta de Andaluc´ıa under grant P12-FQM-1658. References [1] A. Amador, E. Freire, E. Ponce, and J. Ros. On discontinuous piecewise linear models for memristor oscillators. International Journal of Bifurcation and Chaos, 27(06):1730022, jun 2017. [2] H.J. De Blank, Y. Kuznetsov, M.J. Pekk´er, and D.W.M. Veldman. Degenerate Bogdanov- Takens bifurcations in a one-dimensional transport model of a fusion plasma. Physica D: Nonlinear Phenomena, 331:13 – 26, 2016. [3] T.R. Blows and L.M. Perko. Bifurcation of limit cycles from centers and separatrix cycles of planar analytic systems. SIAM Review, 36(3):341376, 1994. 20 A. AMADOR AND E. FREIRE AND E. PONCE [4] H. Chen and X. Chen. Dynamical analysis of a cubic li´enard system with global parameters. Nonlinearity, 28(10):35353562, Apr 2015. [5] H. Chen and X. Chen. Dynamical analysis of a cubic li´enard system with global parameters (ii). Nonlinearity, 29(6):1798, 2016. [6] H. Chen and X. Chen. Global phase portrait of a degenerate bogdanov-takens system with symmetry II. Discrete & Continuous Dynamical Systems-B, 23(10):4141–4170, 2018. [7] H. Chen, X. Chen, and J. Xie. Global phase portrait of a degenerate bogdanov-takens system with symmetry. Discrete & Continuous Dynamical Systems-B, 22(4):1273–1293, 2017. [8] G. Dangelmayr and J. Guckenheimer. On a four parameter family of planar vector fields. Archive for Rational Mechanics and Analysis, 97(4):321352, 1987. [9] A. Dhooge, W. Govaerts, Y. Kuznetsov, H.G.E. Meijer, and B. Sautois. New features of the software matcont for bifurcation analysis of dynamical systems. Mathematical and Computer Modelling of Dynamical Systems, 14(2):147175, 2008. [10] F. Dumortier and C. Li. Perturbations from an elliptic hamiltonian of degree four. Journal of Differential Equations, 175(2):209243, 2001. [11] F. Dumortier and C. Li. Perturbation from an elliptic Hamiltonian of degree four-III global centre. Journal of Differential Equations, 188(2):473511, 2003. [12] F. Dumortier and C. Li. Perturbation from an elliptic Hamiltonian of degree fourIV global centre. Journal of Differential Equations, 188(2):512566, 2003. [13] F. Dumortier, R. Roussarie, and J. Sotomayor. Generic 3-parameter families of vector fields on the plane, unfolding a singularity with nilpotent linear part. the cusp case of codimension 3. Ergodic Theory and Dynamical Systems, 7(03), 1987. [14] Freddy Dumortier. Bifurcations of planar vector fields: nilpotent singularities and Abelian integrals. Springer-Verlag, Berlin New York, 1991. [15] F. Fern´andez, E. Freire, L. Pizarro, and A.J. Rodr´ıguez-Luis. Proc NDES 96: 4th workshop on nonlinear dynamics of electronic systems (seville, spain). ed J.L. Huertas and A.J. Rodr´ıgues- Luis, page 3216, 1996. [16] E. Freire, L. Pizarro, A. J. Rodr´ıguez-Luis, and F. Fern´andez S´anchez. Multiparametric bifurcations in an enzyme-catalyzed reaction model. International Journal of Bifurcation and Chaos, 15(3):905947, 2005. [17] E. Freire, L. Pizarro, and A.J. Rodr´ıguez-Luis. Examples of non-degenerate and degenerate cuspidal loops in planar systems. Dynamics and Stability of Systems, 14(2):129161, 1999. [18] E. Freire, L. Pizarro, and A.J. Rodr´ıguez-Luis. Numerical continuation of homoclinic orbits to non-hyperbolic equilibria in planar systems. Nonlinear Dynamics, 23(4):353375, 2000. [19] P. Godara, D. Dudkowski, A. Prasad, and T. Kapitaniak. New topological tool for multistable dynamical systems. Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(11):111101, nov 2018. [20] J. Guckenheimer, R. Rand, and D. Schlomiuk. Degenerate homoclinic cycles in perturbations of quadratic hamiltonian systems. Nonlinearity, 2(3):405, 1989. [21] Jack Hale. Dynamics and Bifurcations. Springer New York, New York, NY, 1991. [22] M. Itoh and L.O. Chua. Memristor oscillators. International Journal of Bifurcation and Chaos, 18(11):3183–3206, 2008. [23] A.I. Khibnik, B. Krauskopf, and C. Rousseau. Global study of a family of cubic lienard equations. Nonlinearity, 11(6):1505, 1998. [24] L. Kong and C. Zhu. Bogdanov-takens bifurcations of codimensions 2 and 3 in a Leslie-Gower predator-prey model with michaelis-menten-type prey harvesting. Mathematical Methods in the Applied Sciences, 40(18):67156731, Mar 2017. [25] I.A. Korneev and V.V. Semenov. Andronov- Hopf bifurcation with and without parameter in a cubic memristor oscillator with a line of equilibria. Chaos: An Interdisciplinary Journal of Nonlinear Science, 27(8):081104, aug 2017. [26] N.V. Kuznetsov, G.A. Leonov, and S.M. Seledzhi. Hidden oscillations in nonlinear control systems. IFAC Proceedings Volumes, 44(1):2506–2510, jan 2011. [27] G.A. Leonov, N.V. Kuznetsov, and V.I. Vagaitsev. Localization of hidden Chua’s attractors. Physics Letters A, 375(23):2230–2233, jun 2011. [28] V.K. Melnikov. On the stability of a center for time-periodic perturbations. Trans. Moscow Math, 12:1–57, 1963. BIFURCATION SET FOR A DISREGARDED BOGDANOV-TAKENS UNFOLDING 21 [29] M. Messias, C. Nespoli, and V.A. Botta. Hopf bifurcation from lines of equilibria without parameters in memristor oscillators. International Journal of Bifurcation and Chaos, 20(02):437– 450, feb 2010. [30] E. Ponce, A. Amador, and J. Ros. Unravelling the dynamical richness of 3D canonical memristor oscillators. Microelectronic Engineering, 182:15–24, oct 2017. [31] A.J. Rodr´ıguez-Luis, E. Freire, and E. Ponce. A method for homoclinic and heteroclinic continuation in two and three dimensions. Continuation and Bifurcations: Numerical Techniques and Applications, Springer Netherlands, pages 197–210, 1990. [32] S. Schecter. The saddle-node separatrix-loop bifurcation. SIAM Journal on Mathematical Analysis, 18(4):11421156, 1987. [33] V. Varshney, S. Sabarathinam, A. Prasad, and K. Thamilmaran. Infinite number of hidden attractors in memristor-based autonomous duffing oscillator. International Journal of Bifurcation and Chaos, 28(01):1850013, 2018. [34] V. Varshney, S. Sabarathinam, K. Thamilmaran, M.D. Shrimali, and A. Prasad. Existence and control of hidden oscillations in a memristive autonomous duffing oscillator. Nonlinear Dynamical Systems with Self-Excited and Hidden Attractors Studies in Systems, Decision and Control, page 327344, 2018. [35] Stephen Wiggins. Introduction to applied nonlinear dynamical systems and chaos. Springer, 2003. [36] V. Witte, W. Govaerts, Y. Kuznetsov, and M. Friedman. Interactive initialization and continuation of homoclinic and heteroclinic orbits in MATLAB. ACM Transactions on Mathematical Software, 38(3):134, Jan 2012. Appendix: Dimensional reduction in 3D memristor oscillators We consider a family of three-dimensional systems, which is general enough to capture all the mathematical models of memristor oscillators given in (30). Such family was been studied in [1] and [30], where the authors showed that the dynamics of such a family of three-dimensional systems is essentially ruled by a one parameter set of two-dimensional systems. We consider the system (46) ˙x=a11W(z)x+a12y, ˙y=a21x+a22y, ˙z=x, where the constants a11, a12, a21, a22 ∈Rand the function Wallows to define a continuous function (47) q(z) = Zz 0 W(s)ds. The next result guarantees that the dynamics of system (46) is essentially twodimensional, see [1] for a proof. Proposition 16. Consider system (46) where the functions Wand qare related as in (47). For any h∈R, the set (48) Sh={(x, y, z)∈R3:−a22x+a12y−a12a21z+a11a22q(z) = h} is an invariant manifold for the system. Therefore, the system has an infinite family of invariant manifolds foliating the whole R3, and so the dynamics is essentially two-dimensional. In the following result we show that on each invariant set Shgiven in (48), and for any continuous function qdefined as in (47), the dynamics is topologically equivalent to a Li´enard system. Furthermore, we give for any solution of the Li´enard system with a given value of h, the corresponding solution of the 3D canonical model (46). 22 A. AMADOR AND E. FREIRE AND E. PONCE This result is a generalization of Theorem 3 given in [1], where the function qwas considered to be a continuous piecewise linear function. Proposition 17. Consider system (46) with the function qdefined as in (47). If a12 6= 0, then on each invariant set Shgiven in (48), the dynamics is topologically equivalent to the Li´enard system (49) ˙ X=Y−F(X),˙ Y=−g(X) + h, where Fand gare given by (50) F(X) = −a11q(X)−a22X, g(X) = a11a22q(X)−a12a21X Moreover, (X(τ), Y (τ)) ∈R2is a solution of the Li´enard system (49) for a given h∈R, if and only if Eh(X(τ), Y (τ)) ∈R3is a solution of system (46) on Sh, where (51) Eh(X(τ), Y (τ)) =   Y(τ)−F(X(τ)) 1 a12 [(a2 22 +a12a21)Y(τ)−a22Y(τ) + h] X(τ)  . Proof. First, with a12 6= 0 the change of variables (52) x=x, y =a22x−a12y, z =z transforms system (46) into the system ˙ x=f1(z)x−y,(53) ˙ y=f2(z)x, ˙ z=x, where the functions f1and f2are defined as (54) f1(z) = a11W(z) + a22, f2(z) = a22a11W(z)−a12a21. From Proposition 16, the invariant manifolds (48) for system (53)-(54) can be written in the new variables as (55) e Sh={(x, y, z)∈R3:−y+g(z) = h}. Now, replacing the condition given in (55) in the first equation of (53) and removing the unnecessary second equation, we obtain the system (56) ˙ x=f1(z)x−g(z) + h, ˙ z=x. where the function gis defined by (57) g(u) = a11a22q(u)−a12a21u. After the change of variables (58) X=z, Y=−˜ F(z) + x, where Fis (59) ˜ F(z) = a11q(z) + a22z, we obtain ˙ X=˙ z=x=Y+˜ F(X) = Y−(−˜ F(X)), BIFURCATION SET FOR A DISREGARDED BOGDANOV-TAKENS UNFOLDING 23 so that ˙ Y=−˜ F′(z)˙ z+˙ x=−(a11q′(z) + a22) + (f1(z)x−g(z) + h) = =−f1(z)x+f1(z)x−g(z) + h=−g(z) + h, and taking F(X) = −˜ F(X) we obtain system (49)-(50). If (X(τ), Y (τ)) ∈R2is a solution of system (49)-(50) for a given h∈R, we have from (58) that x(τ) z(τ)=Y(τ)−F(z(τ)) X(τ) is a solution of system (56). From (55), we obtain on e Shthat y=g(z)−h, with g as in (57). Thus,  x(τ) y(τ) z(τ)  = Y(τ)−F(X(τ)) g(X(τ)) −h X(τ)  , is a solution of system (53) on Sh. Finally, from (52) we obtain for system (46) the solution x(τ) = x(τ), y(τ) = 1 a12 [a22x(τ)−y(τ)] = 1 a12 [a22Y(τ)−a22F(X(τ)) −g(X(τ)) + h] =1 a12 ha22Y(τ) + a22 ˜ F(X(τ)) −g(X(τ)) + hi, and z(τ) = z(τ). The conclusion follows from the fact that for all Xwe have a22 ˜ F(X)−g(X) = (a2 22 +a12a21)X.  In order to apply the analysis performed to system (3), in what follows we consider the function qdefined by a cubic polynomial, that is, we assume (60) W(z) = 3cz2+ 2az +b, q(z) = cz3+az2+bz, with c6= 0. As a direct consequence of Propositions 16 and 17, we obtain the next result. Corollary 18. Consider system (46) with the functions qand Wdefined as in (60). If a12 6= 0, then on each invariant set Shgiven by Sh={(x, y, z)∈R3:−a22x+a12y+a11a22cz3+aa11a22z2+(ba11a22 −a12a21)z=h} the dynamics is topologically equivalent to the Li´enard system ˙x=y+ca11x3+aa11x2+ (ba11 +a22)x, ˙y=−a11a22cx3−a11a22ax2+ (a12a21 −a11a22b)x+h. (61) Moreover, (x(τ), y (τ)) ∈R2is a solution of the Li´enard system (61) for a given h∈R, if and only if Eh(x(τ), y (τ)) ∈R3is a solution of system (46) on Sh,where (62) Eh(x(τ), y (τ)) =   y(τ) + ca3 11x(τ)3+aa2 11x(τ)2+ (ba11 +a22)x(τ)2 1 a12 [(a2 22 +a12a21)y(τ)−a22y(τ) + h] x(τ)  . In the next Proposition, we show that system (61) can be written into the form (1). 24 A. AMADOR AND E. FREIRE AND E. PONCE Proposition 19. The following statements hold for system (61). (a) If a22 6= 0 and a11a22 <0then the system can be written into the form (63) ˙x=y, ˙y=µ1+µ2x+cx3+µ3y+ 3ca11x2y. where the new parameters µ1, µ2and µ3are given by µ1=27ch +a11a22a(9cb −2a2)−9caa12a21 27c2(−a11a22)5/2, µ2=a11a22(a2−3cb) + 3ca12a21 3c(a11a22)2, µ3=a11(a2−3cb)−3ca22 3ca11a22 . (64) (b) If a22 = 0 then the system can be written into the form (65) ˙x=y, ˙y=µ1+µ2x+µ3y+ 3ca11x2y, where the new parameters µ1, µ2and µ3are defined by (66) µ1=h−aa12a21 3c, µ2=a12a21, µ3=ba11 −a2a11 3c. Proof. First, the change of variables u=x+a 3c, v =y+2 27 a3 c2a11 −1 3 a ca22 −1 3ab ca11, transforms system (61) into ˙u=v+ca11u3+λ1u, ˙v=−ca11a22u3+λ2u+λ3, (67) where the new parameters are λ1=a22 +ba11 −1 3 a2 ca11, λ2=a12a21 −ba11a22 +1 3 a2 ca11a22, λ3=h+1 3ab ca11a22 −1 3 a ca12a21 −2 27 a3 c2a11a22. (68) If a11a22 <0,the change of variable x=1 (−a11a22)1/2u, y =v, τ =1 −a11a22 t, transforms system (67)-(68) into ˙x=1 (−a11a22)3/2y+ca11x3−λ1 a11a22 x, ˙y=λ3 −a11a22 +λ2 (−a11a22)1/2x+c(−a11a22)3/2x3 and taking into account that ¨x=1 (−a11a22)3/2˙y+ 3ca11x2˙x−λ1 a11a22 ˙x, and after some algebra, statement (a) follows. If a22 = 0,then from system (67), we obtain statement (b) after a direct computation. 