# Bifurcation Analysis of a Van der Pol-Duffing Circuit with Parallel Resistor

Authors: Denis de Carvalho Braga; Luis Fernando Mello; Marcelo Messias
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

We study the local codimension one, two, and three bifurcations which
occur in a recently proposed Van der Pol-Duffing circuit (ADVP) with parallel
resistor, which is an extension of the classical Chua's circuit with cubic nonlinearity.
The ADVP system presents a very rich dynamical behavior, ranging from
stable equilibrium points to periodic and even chaotic oscillations. Aiming to contribute
to the understand of the complex dynamics of this new system we present
an analytical study of its local bifurcations and give the corresponding bifurcation
diagrams. A complete description of the regions in the parameter space for
which multiple small periodic solutions arise through the Hopf bifurcations at the
equilibria is given. Then, by studying the continuation of such periodic orbits,
we numerically find a sequence of period doubling and symmetric homoclinic bifurcations
which leads to the creation of strange attractors, for a given set of the
parameter values.


---

## Body

## 1. Introduction and Statement of the Main Results

In this paper we study the local codimension one, two, and three bifurcations and the respective qualitative  changes in the dynamics of the following system of nonlinear equations:(1.1)ẋ=dxdτ=-ν(x3-μx-y),ẏ=dydτ=x-αy-z,ż=dzdτ=βy,
where (x,y,z)∈ℝ3 are the state variables and μ∈ℝ, α≥1, ν>0, and β>0 are real parameters.As far as we  know, system (1.1) was proposed and firstly studied in [1], and it can be obtained from the system(1.2)c1v1′=-[iN+(v1-v2)/R],c2v2′=(v1-αv2)/R-iL,iL′=v2/L,
by the following changes in variables, parameters, and rescaling in time:(1.3)x=bRv1,y=bRv2,z=bR3iL,μ=-(1+aR),ν=c2/c1,β=c2R2/L,τ=t/Rc2.In (1.2) v1 and v2 are the voltages across the capacitors C1 and C2, respectively, iN is the current through the nonlinear diode N, iL is the current through the inductor, c1 and c2 are the capacitances, L is the inductance, R is the linear resistor, Rp is the parallel resistor, α=1+(R/Rp), and the prime denotes derivatives with respect to the time t. It is assumed that the current through the nonlinear diode is given by iN=iN(v1)=av1+bv13, with a<0 and b>0. See Figure 1 and [1] for more details.Figure 1
ADVP circuit with parallel resistorRp.Despite the simplicity of the electronic circuit shown in Figure1, the related system (1.1) has a rich dynamical behavior, ranging from stable equilibrium points to periodic and even chaotic oscillations, depending on the parameter values. The study of simple nonlinear electronic circuits presenting chaotic behavior was initiated  by  Chua in the mid-1980s (see [2]) and the interest on the subject has grown in the last decades. In fact, it has useful applications like chaos synchronization which is applied in industry and secure communications, beyond other areas of sciences (see [1] and references therein and also [3, 4]). System (1.1) with α=1 reduces to a system equivalent to the classical Chua's differential equations with cubic nonlinearity. The authors have recently developed a complete study of degenerate Hopf bifurcations for this case in [5], answering a  challenge proposed by  Moiola and Chua  in [6].System (1.1) has the equilibrium point E0=(0,0,0), which exists for any parameter values. For μ>0 it also has the symmetric equilibria E±=(±μ,0,±μ). Codimension one Hopf bifurcation which occurs at the equilibrium point E0 was studied in [1], by using a result of Hsü and Kazarinoff in [7].In this paper by using the classical projection method which allows one the calculation of the Lyapunov coefficients associated to the Hopf bifurcations we study all possible bifurcations (generic and degenerate ones) which occur at the equilibriaE0 and  E± of system (1.1). In this way the analyses presented in [1] are extended and completed. More precisely, we prove the  following statements.(a)
For the equilibriumE0=(0,0,0) the Hopf hypersurface H1 is obtained in the space of parameters (μ,ν,α,β) and the first Lyapunov coefficient l1 is calculated. It is shown that this coefficient vanishes on a 2-dimensional surface contained in H1, giving rise to codimension two Hopf bifurcations (see Figure 2(a)). Then the second Lyapunov coefficient l2 is calculated and it is established that this coefficient is always negative on the surface  {l1=0}.(b)
For the symmetric equilibriaE±=(±μ,0,±μ) the Hopf hypersurface H2 is obtained in the space of parameters (μ,ν,α,β), the first Lyapunov coefficient l1 is calculated, and it is shown that this coefficient vanishes on a 2-dimensional surface contained in H2 (see Figure 2(b)), giving rise to codimension two bifurcations. The second Lyapunov coefficient l2 is calculated and it is established that this coefficient also vanishes on a 2-dimensional surface contained in H2 (see Figure 2(b)), giving rise to codimension three bifurcations. The third Lyapunov coefficients along the curve C given by the intersection of the surfaces {l1=0} and {l2=0} are obtained and found to be positive.Surfaces{l1=0} for E0 (a) and E± (b), C={l1=0}∩{l2=0}.
(a)(b)The corresponding bifurcation diagrams are traced for the above bifurcations (see Figures3, 4, 7, and 8 of Section 4). From statement (a) it follows that the maximum number of small periodic orbits bifurcating from the origin is 2. Furthermore, from statement (b) one can deduce that there is a region in the parameter space for which an attracting periodic orbit and two other unstable periodic orbits coexist with the attracting equilibrium points E±. See Figures 7 and 8 in Section 4.(a) Typical pointsN1, N2, N3, N4, N5, and N6 in the parameter space near the point R.  (b) Phase portraits of system (1.1) for the flow restricted to the center manifold and its continuation. For parameters at N1 the equilibrium is asymptotically stable; for parameters at N2 the equilibrium is a  weak stable focus (Hopf point with l1 negative); for parameters at N3 the equilibrium is unstable and a stable limit cycle appears from a Hopf bifurcation; for parameters at N4 the equilibrium is an weak unstable focus (Hopf point with l1 positive) and there is a stable limit cycle; for parameters at N5 the equilibrium is asymptotically stable and an unstable limit cycle appears from a Hopf bifurcation, so there are two limit cycles encircling the equilibrium; for parameters at N6 the equilibrium is asymptotically stable and the two cycles collide, giving rise to a nondegenerate fold bifurcation of the cycles.
(a)(b)(a) Typical pointsM1, M2, M3, M4, M5, and M6 in the parameter space near the point T.  (b) Phase portraits of system (1.1) for the flow restricted to the center manifold and its continuation. For parameters at M1 the equilibrium is unstable; for parameters at M2 the equilibrium is an weak unstable focus (Hopf point with l1 positive); for parameters at M3 the equilibrium is stable and an unstable limit cycle appears from a Hopf bifurcation; for parameters at M4 the equilibrium is an weak stable focus (Hopf point with l1 negative) and there is an unstable limit cycle; for parameters at M5 the equilibrium is unstable and a stable limit cycle appears from a Hopf bifurcation, so there are two limit cycles encircling the equilibrium; for parameters at M6 the equilibrium is unstable and the two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles.
(a)(b)We go further and perform a numerical study of the continuation of periodic orbits which arise from the Hopf bifurcations at the pointsE±. Through this analysis we have found period doubling and homoclinic bifurcations which seems to lead to the creation of strange attractors of system (1.1), for some parameter values (see Section 5).The rest of this paper is organized as follows. In Section2 through a linear analysis of system (1.1) we obtain the Hopf hypersurfaces for E0 and E±. In Section 3 following [8–10] we present a brief review of the methods used to study codimension one, two, and three Hopf bifurcations, describing in particular how to calculate the Lyapunov coefficients related to the stability of the equilibrium point as well as of the periodic orbits which appear in these bifurcations. In general the Lyapunov coefficients are very difficult to be obtained. These methods are used in Section 4 to prove the main results of this paper, described in statements (a) and (b). In Section 5 we present numerical simulations for particular values of the parameters which illustrate and corroborate the analytical results obtained. Finally, in Section 6 we make some concluding remarks.


## 2. Linear Analysis of System (1.1)

In a vectorial notation which will be useful in the next calculations system (1.1) can be written as(2.1)x′=f(x,ζ)=(-ν(x3-μx-y),x-αy-z,βy),
where(2.2)x=(x,y,z)∈ℝ3,ζ=(μ,ν,α,β)∈(-∞,∞)×(0,∞)×[1,∞)×(0,∞).The originE0 is an equilibrium point of system (1.1) for all values of the parameters. The characteristic polynomial of the Jacobian matrix of the function f given in (2.1) at E0 is given by(2.3)p(λ)=λ3+(α-μν)λ2+(β-ν-αμν)λ-μβν.
If μ>0 then the term -μβν in the above polynomial is negative and this implies that E0 is an unstable equilibrium point.For the sake of completeness we state the following lemma (Routh-Hurwitz stability criterion) whose proof can be found in [11, page 58].Lemma 2.1.
The polynomialL(λ)=p0λ3+p1λ2+p2λ+p3,  p0>0, with real coefficients, has all roots with negative real parts if and only if the numbers p1, p2, p3 are positive and the inequality p1p2>p0p3 is satisfied.The following proposition is a direct consequence of Lemma2.1.Proposition 2.2.
Considerμ<0 in system (1.1). If αμ≤-1 then the equilibrium point  E0 is asymptotically stable for all β>0 and ν>0. If αμ>-1 and
(2.4)β>β1c=ν(α-μν)(1+αμ)α,
then system (1.1) has an asymptotically stable equilibrium point at E0. If αμ>-1 and 0<β<β1c then E0 is unstable.Forμ=0 the origin becomes a nonhyperbolic equilibrium point and the situation is highly degenerate. This case, which will not be considered in this work, can be studied by using the methods presented for instance in [12].Forμ>0 two new equilibria E±=(±μ,0,±μ) appear. System (1.1) is invariant under the change of coordinates (x,y,z)→(-x,-y,-z). So the stability of the equilibrium E- can be obtained from the stability of E+.The characteristic polynomial of the Jacobian matrix of the functionf given in (2.1) at E+ is given by(2.5)p(λ)=λ3+(α+2μν)λ2+(β-ν+2αμν)λ+2μβν.
The following proposition is also a direct consequence of Lemma 2.1.Proposition 2.3.
Considerμ>0 in system (1.1). If  2αμ≥1 then the equilibrium point E+ is asymptotically stable for all β>0 and ν>0. If 2αμ<1 and
(2.6)β>β2c=ν(α+2μν)(1-2αμ)α,
then system (1.1) has an asymptotically stable equilibrium point at E+. If 2αμ<1 and 0<β<β2c then E+ is unstable.The equationsβ=β1c and β=β2c in (2.4) and (2.6) give the equations of the Hopf hypersurfaces H1 and H2 in the parameter space (μ,ν,α,β). These equations will be used in Section 4 in the study of degenerate Hopf bifurcations which occur at the equilibria E0 and E± of system (1.1).


## 3. Outline of the Hopf  Bifurcation Methods

This section is a review of the projection method described in [8] for the calculation of the first and second Lyapunov coefficients associated to Hopf bifurcations, denoted by  l1 and l2, respectively. This method was extended to the calculation of the third and fourth Lyapunov coefficients in [9, 10]. Other equivalent definitions and algorithmic procedures to write the expressions of the Lyapunov coefficients for two-dimensional systems can be found in Andronov et al. [13] and Gasull and Torregrosa [14], among others, and can be adapted to the three-dimensional system of this paper if it is restricted to the center manifold. See also [15].Consider the differential equation(3.1)x′=f(x,ζ),
where x∈ℝ3 and ζ∈ℝ4 are, respectively, vectors representing phase variables and control parameters. Assume that f is of class C∞ in ℝ3×ℝ4. Suppose that (3.1) has an equilibrium point x=x0 at ζ=ζ0 and, denoting the variable x-x0 also by x, write(3.2)F(x)=f(x,ζ0)
as(3.3)F(x)=Ax+12B(x,x)+16C(x,x,x)+124D(x,x,x,x)+1120E(x,x,x,x,x)+1720K(x,x,x,x,x,x)+15040L(x,x,x,x,x,x,x)+O(∥x∥8),
where A=fx(0,ζ0) and, for i=1,2,3,(3.4)Bi(x,y)=∑j,k=13∂2Fi(ξ)∂ξj∂ξk|ξ=0xjyk,Ci(x,y,z)=∑j,k,l=13∂3Fi(ξ)∂ξj∂ξk∂ξl|ξ=0xjykzl,
and so on for Di, Ei, Ki, and Li.Suppose that(x0,ζ0) is an equilibrium point of (3.1) where the Jacobian matrix A has a pair of purely imaginary eigenvalues λ2,3=±iω0, ω0>0, and admits no other eigenvalue with zero real part. Let Tc be the generalized eigenspace of A corresponding to λ2,3. By this it is meant the largest subspace invariant by A on which the eigenvalues are λ2,3.Letp,q∈ℂ3 be vectors such that(3.5)Aq=iω0q,A⊤p=-iω0p,〈p,q〉=∑i=13p̅iqi=1,
where A⊤ is the transpose of the matrix A. Any vector y∈Tc can be represented as y=wq+w̅q̅, where w〈p,y〉∈ℂ. The two-dimensional center manifold associated to the eigenvalues λ2,3=±iω0 can be parameterized by the variables w and w̅ by means of an immersion of the form x=H(w,w̅), where H:ℂ2→ℝ3 has a Taylor expansion of the form(3.6)H(w,w̅)=wq+w̅q̅+∑2≤j+k≤71j!k!hjkwjw̅k+O(|w|8),
with hjk∈ℂ3 and hjk=h̅kj. Substituting this expression into (3.1) we obtain the following differential equation:(3.7)Hww′+Hw̅w̅′=F(H(w,w̅)),
where F is given by (3.2). The complex vectors hij are obtained solving the system of linear equations defined by the coefficients of (3.7), taking into account the coefficients of F (see   [9, Remark  3.1, page 27]), so that system (3.7), on the chart w for a central manifold, writes as follows:(3.8)w′=iω0w+12G21w|w|2+112G32w|w|4+1144G43w|w|6+O(|w|8),
with Gjk∈ℂ.Thefirst Lyapunov coefficientl1 is defined by(3.9)l1=12ReG21,
where G21=〈p,ℋ21〉, and ℋ21=C(q,q,q̅)+B(q̅,h20)+2B(q,h11).Thesecond Lyapunov coefficient is defined by(3.10)l2=112ReG32,
where G32=〈p,ℋ32〉. The expression for ℋ32 can be found in   [9, equation (36), page 28].Thethird Lyapunov coefficient is defined by(3.11)l3=1144ReG43,
where G43=〈p,ℋ43〉. The expression for ℋ43 can be found in   [9, equation (44), page 30].AHopf point(x0,ζ0) of system (3.1) is an equilibrium point where the Jacobian matrix A has a pair of purely imaginary eigenvalues λ2,3=±iω0,  ω0>0, and the other eigenvalue λ1≠0. From the Center Manifold Theorem, at a Hopf point a two-dimensional center manifold is well defined, it is invariant under the flow generated by (3.1) and can be continued with arbitrary high class of differentiability to nearby parameter values (see [8]).A Hopf point is calledtransversal if the parameter dependent complex eigenvalues cross the imaginary axis with nonzero derivative. In a neighborhood of a transversal Hopf point with l1≠0 the dynamic behavior of the system (3.1), reduced to the family of parameter-dependent continuations of the center manifold, is orbitally topologically equivalent to the following complex normal form w′=(η+iω)w+l1w|w|2, where w∈ℂ, η, ω, and l1 are real functions having derivatives of arbitrary higher order, which are continuations of 0, ω0, and the first Lyapunov coefficient at the Hopf point. See [8, 9] for details. When l1<0(l1>0) one family of stable (unstable) periodic orbits can be found on this family of manifolds, shrinking to an equilibrium point at the Hopf point.AHopf point of codimension 2 is a Hopf point where l1 vanishes. It is called transversal if {η=0} and {l1=0} have transversal intersections, where η=η(ζ) is the real part of the critical eigenvalues. In a neighborhood of a transversal Hopf point of codimension 2 with l2≠0 the dynamic behavior of the system (3.1)  reduced to the family of parameter-dependent continuations of the center manifold is orbitally topologically equivalent to w′=(η+iω0)w+τw|w|2+l2w|w|4, where η and τ are unfolding parameters. See [8, 9]. The bifurcation diagrams for l2≠0 can be found in [8, page 313] and in [16].AHopf point of codimension 3 is a Hopf point of codimension 2 where l2 vanishes. It is called transversal if {η=0}, {l1=0}, and {l2=0} have transversal intersections. In a neighborhood of a transversal Hopf point of codimension 3 with l3≠0 the dynamic behavior of the system (3.1), reduced to the family of parameter-dependent continuations of the center manifold, is orbitally topologically equivalent to w′=(η+iω0)w+τw|w|2+νw|w|4+l3w|w|6, where η, τ,  and ν are unfolding parameters. The bifurcation diagram for l3≠0 can be found in Takens [16].


## 4. Hopf Bifurcations in System (1.1)

### 4.1. Hopf Bifurcations atE0

In this subsection we study the stability ofE0 under the conditions β=β1c and -1<αμ<0, that is, on the Hopf hypersurface H1 complementary to the range of validity of Proposition 2.2.The Jacobian matrix of the functionf given in (2.1) at E0 is given by(4.1)A=Df(E0)=(μνν01-α-10β1c0).
Then, using the notation of the previous section (see expression (3.3)) the multilinear symmetric functions corresponding to f can be written as(4.2)B(x,y)=D(x,y,z,u)=E(x,y,z,u,v)=(0,0,0),C(x,y,z)=(-6νx1y1z1,0,0).The eigenvalues ofA are(4.3)λ1=-α+μν<0,λ2,3=±iω0=±iν-μ(1+αμ)α,
and from (3.5) one has
(4.4)q=(-iω0νμν-iω0,iω0,β1c),p̅=1ρ0(iω0(μν-iω0),-iω0(μν-iω0)2,(μν-iω0)2),
where
(4.5)ρ0=β1c(μν-iω0)2+(ν+(μν-iω0)2)ω02.Using these calculations we prove the next two theorems, from which statement (a) of the introduction follows.Theorem 4.1.
Consider the four-parameter family of differential equations (1.1). The first Lyapunov coefficient associated to the equilibrium E0 is given by
(4.6)l1(μ,ν,α,β1c)=-3α(αμ+1)ν3(α2-2μνα-ν)2(α3-2μνα2-μν2).
If
(4.7)s(μ,ν,α)=-α2+2μνα+ν
is different from zero then system (1.1) has a transversal Hopf point at E0 for β=β1c and -1<αμ<0.
More precisely, if(μ,ν,α,β1c)∈S∪U (see Figure 2(a)), then system (1.1) has a Hopf point of codimension 1 at E0. If (μ,ν,α,β1c)∈S then E0 is asymptotically stable and for each β<β1c, but close to β1c, there exists a stable periodic orbit near the unstable equilibrium point E0. If (μ,ν,α,β1c)∈U then E0 is unstable and for each β>β1c, but close to β1c, there exists an unstable periodic orbit near the asymptotically stable equilibrium point E0.Proof.
As the functionB(x,y)=(0,0,0) then the expression G21 in (3.9) reduces to G21〈p,C(q,q,q̅)〉. Performing the calculations we get
(4.8)G21=-6ν4ω04(μ2ν2+ω02)(β1c(μν-iω0)2+ω02(μ2ν2-2iμω0ν+ν-ω02)).
Substituting the expressions of ω0 and β1c into the above expression of G21 and taking the real part we have, from (3.9), the expression of the first Lyapunov coefficient l1 given in  (4.6). Note that the sign of l1 is given by the sign of the function s defined in (4.7).
It remains only to verify the transversality condition of the Hopf bifurcation. In order to do so, consider the family of differential equations (1.1) regarded as dependent on the parameter β. The real part, η=η(β), of the pair of complex eigenvalues at the critical parameter β=β1c verifies
(4.9)η′(β1c)=Re〈p,dAdβ|β=β1cq〉=-α22(α3-2μνα2-μν2)<0.
Therefore, the transversality condition at the Hopf point holds.
The setsS and U illustrated in Figure 2(a) are defined implicitly as {l1<0} and {l1>0},  respectively. The theorem is proved.Theorem 4.2.
Consider the four-parameter family of differential equations (1.1) restricted to β=β1c  and -1<αμ<0. Then the second Lyapunov coefficient l2 associated to the equilibrium E0 is negative on the Hopf hypersurface where the first Lyapunov coefficient l1 vanishes.Proof.
Following the notation introduced in Section3, we have B(x,y)=D(x,y,z,u)=E(x,y,z,u,v)=(0,0,0). Then
(4.10)ℋ32=3C(q,q,h̅21)+6C(q,q̅,h21)+3C(q̅,q̅,h30).
The expressions of the complex vectors h21 and h30 are too long to be put in print. The interested author can encounter such expressions and all the details of the calculations presented in this paper in the website [17].By direct calculations the second Lyapunov coefficientl2 associated to E0 can be written as(4.11)l2(μ,ν,α,β1c)=-3α2(αμ+1)2ν5N(μ,ν,α,β1c)2μ(-α3+2μνα2+μν2)3(-α3+2μνα2+8μ2ν2α+9μν2),
where N(μ,ν,α,β1c) has the form(4.12)α10-8μνα9+ν(16μ2ν-1)α8+2μν2(8μ2ν-7)α7+2μ2ν3(61-40μ2ν)α6+2μν3(32ν2μ4-112νμ2+13)α5-μ2ν4(64νμ2+19)α4+40μ3ν5(10μ2ν-9)α3+μ2ν5(-192ν2μ4+724νμ2-121)α2-8μ3ν6(43μ2ν-34)α-142μ4ν7.
Along the surface {l1=0} we have, from   (4.7) of Theorem 4.1, s(μ,ν,α)=-α2+2μνα+ν=0, or equivalently, να2/(1+2αμ). Substituting this expression into (4.11) one has(4.13)l2|{l1=0}=-9α11(αμ+1)(4αμ-1)(2αμ+1)4(8αμ-1)<0.Here the transversality condition is equivalent to the transversality of the hypersurfacesβ=β1c and  l1(μ,ν,α,β)=0. The bifurcation diagram for a typical point R on the surface {l1=0} is depicted in Figure 3.


### 4.2. Hopf Bifurcations atE±

In this subsection we study the stability ofE± under the conditions β=β2c and 0<2αμ<1, that is, on the Hopf hypersurface H2 complementary to the range of validity of Proposition 2.3.The Jacobian matrix of the functionf given in (2.1) calculated at E+ is given by(4.14)A=Df(E+)=(-2μνν01-α-10β2c0).
Then, using the notation of Section 3 (see expression (3.3)) the multilinear symmetric functions corresponding to f can be written as(4.15)B(x,y)=(-6νμx1y1,0,0),C(x,y,z)=(-6νx1y1z1,0,0),D=E=K=L=(0,0,0).The eigenvalues ofA are(4.16)λ1=-α-2μν<0,λ2,3=±iω0=±iν2μ(1-2αμ)α,
and from (3.5) one has(4.17)q=(iω0ν2μν+iω0,iω0,β2c),(4.18)p̅=1ρ0(-iω0(2μν+iω0),-iω0(2μν+iω0)2,(2μν+iω0)2),
where(4.19)ρ0=β2c(2μν+iω0)2+(ν+(2μν+iω0)2)ω02.Using these calculations we prove the next three theorems, from which statement (b) of introduction follows.Theorem 4.3.
Consider the four-parameter family of differential equations (1.1). The first Lyapunov coefficient associated to the equilibrium E+ is given by
(4.20)l1(μ,ν,α,β2c)=-3α(2αμ-1)ν3g(μ,ν,α)(α3+4μνα2+2μν2)(α3+4μνα2-12μ2ν2α+8μν2),
where
(4.21)g(μ,ν,α)=α5+8μνα4+ν(4μ2ν-1)α3+2μν2(5-24μ2ν)α2+80μ2ν3α-20μν3.
If g(μ,ν,α) is different from zero then system (1.1) has a transversal Hopf point at E+ for β=β2c and 0<2αμ<1.
More precisely, if(μ,ν,α,β2c)∈S∪U (see Figure 2(b)) then system (1.1) has a Hopf point of codimension 1 at E+. That is, if (μ,ν,α,β2c)∈S then E+ is asymptotically stable and for each β<β2c, but close to β2c, there exists a stable periodic orbit near the unstable equilibrium point E+; if (μ,ν,α,β2c)∈U then E+ is unstable and for each β>β2c, but close to β2c, there exists an unstable periodic orbit near the asymptotically stable equilibrium point E+.Proof.
The theorem follows by expanding the expressions in definition of the first Lyapunov coefficient (3.9). It relies on extensive calculations involving the vector q in (4.17), the vector p in (4.18) and the multilinear functions B and C. The inclusion here of the enormous expressions obtained from the mentioned calculations would not bring any contribution in clarifying the reading of the text. Then the detailed calculations related to this proof, corroborated by Computer Algebra, have been posted in [17], including its implementation using the software MATHEMATICA 5 [18].
Notice that the sign of the first Lyapunov coefficient is given by the sign of the function-g defined in (4.21) since the denominator of l1 and the expression 3(2αμ-1)ν3 are negative. In Figure 2(b) the surface {l1=0} and the regions S and U are illustrated. These regions are defined implicitly as {l1<0} and {l1>0}, respectively.
It remains only to verify the transversality condition of the Hopf bifurcation. In order to do so, consider the family of differential equations (1.1) regarded as dependent on the parameter β. The real part, η=η(β), of the pair of complex eigenvalues at the critical parameter β=β2c verifies
(4.22)η′(β2c)=Re〈p,dAdβ|β=β2cq〉=-α22(α3+4μνα2+2μν2)<0.
Therefore, the transversality condition holds at the Hopf point. The theorem is proved.Theorem 4.4.
Consider the four-parameter family of differential equations (1.1) restricted to β=β2c and 0<2αμ<1. Then the second Lyapunov coefficient l2 associated to the equilibrium E+ is given by
(4.23)l2(μ,ν,α,β2c)=α2(2αμ-1)ν5h(μ,ν,α)4μ(α3+4μνα2-32μ2ν2α+18μν2)d(μ,ν,α),
where the function h(μ,ν,α) is given by
(4.24)21μα20+(588μ2ν-4)α19+3μν(1876μ2ν-53)α18+4ν(3192ν2μ4-294νμ2+1)α17+μν2(-121296ν2μ4+14156νμ2-65)α16-12μ2ν3(59248ν2μ4-18704νμ2+473)α15+μν3(114240ν3μ6+716976ν2μ4-62932νμ2+199)α14+16μ2ν4(477120ν3μ6-165272ν2μ4+769νμ2+119)α13+4μ3ν5(1634304ν3μ6-4121136ν2μ4+751264νμ2-16945)α12-16μ2ν5(2225664ν4μ8-327552ν3μ6-585624ν2μ4+53587νμ2-262)α11-4μ3ν6(6967296ν4μ8-29565184ν3μ6+7648944ν2μ4+81592νμ2-11129)α10+16μ4ν7(4644864ν4μ8+773376ν3μ6-8268704ν2μ4+1627004νμ2-41405)α9+4μ3ν7(-80474112ν4μ8+32219456ν3μ6+13655344ν2μ4-1903712νμ2+17531)α8+16μ4ν8(5308416ν4μ8+34114176ν3μ6-13584144ν2μ4-51152νμ2+39713)α7+32μ5ν9(1327104ν4μ8-10505088ν3μ6-13899000ν2μ4+4382138νμ2-130059)α6-16μ4ν9(11446272ν4μ8-35076672ν3μ6-10290192ν2μ4+2504248νμ2-37769)α5+32μ5ν10(10175616ν3μ6-15588272ν2μ4-317442νμ2+129593)α4-64μ6ν11(4720704ν2μ4-3846404νμ2+142349)α3+128μ5ν11(1199432ν2μ4-495365νμ2+13321)α2-256μ6ν12(157944μ2ν-25841)α+4308992μ7ν13,
and the function d(μ,ν,α) in the denominator is given by
(4.25)(α6+8μνα5+4μ2ν2α4+2μν2(5-24μ2ν)α3+40μ2ν3α2-24μ3ν4α+16μ2ν4)3.
If (μ,ν,α,β2c)∈C1∪C2 (see Figure 2(b)), then system (1.1) has a transversal Hopf point of codimension 2 at E+. More specifically, if (μ,ν,α,β2c)∈C1 then the Hopf point at E+ is asymptotically stable and the bifurcation diagram is illustrated in Figure 3 for a typical point R. If (μ,ν,α,β2c)∈C2 then the Hopf point at E+ is unstable and the bifurcation diagram is drawn in Figure 4 for a typical point T.Proof.
The theorem follows by expanding the expressions in definition (3.10) of the second Lyapunov coefficient. It relies on extensive calculation involving the vector q in (4.17), the vector p in (4.18) and the multilinear functions B, C, D, and E. The calculations in this proof, corroborated by Computer Algebra, have been posted in [17].
In Figure2 we present a geometric synthesis interpreting the long calculations involved in this proof. The sign of -h(μ,ν,α) gives the sign of the second Lyapunov coefficient since the function d(α,μ,ν) in the denominator of l2 is positive. The graph of h(μ,ν,α)=0 on the surface {l1=0}, which determines the curve C, and the signs of the first and second Lyapunov coefficients are also illustrated. As follows, l2<0 on the open region on the surface {l1=0}, denoted by C1. On this region a typical reference point R is depicted. Also l2>0 on the open region on the surface {l1=0}, denoted by C2. This region contains the typical reference point, denoted by T.
The bifurcation diagrams of system (1.1) at the points R and T are illustrated in Figures 3 and 4, respectively. The theorem is proved.Remark 4.5.
There are numerical evidences that the sign of the third Lyapunov coefficient is always positive along the curveC given by the intersection of the surfaces {l1=0} and {l2=0} (see Figure 2(b)). For α=1 see the article in [5]. For α=2 see Theorem 4.6 in what follows. Figure 5 shows the behavior of the third Lyapunov coefficient as a function of α, for (μ,ν,β) on the curve C. So the bifurcation diagrams depicted in Figures 7 and 8 are essentially the same for all values of α≥1.Figure 5
Graph ofl3 with respect to α.All the calculations presented here are illustrated in Figure6 for the cases α=2 (a) and α=3 (b). In this figure are depicted the surfaces β=β2c, the solid curves where l1=0, the regions S and U where l1<0 and l1>0, respectively, the dashed curves where l2=0 and the arcs C1 and C2 where l2<0 and l2>0, respectively. The bifurcation diagrams for typical points R1 and R2 are depicted in Figure 3 while the bifurcation diagrams for typical points T1 and T2 are shown in Figure 4. In the next theorem we analyze the sign of the third Lyapunov coefficient at the point Q1, that is, for the case α=2 or equivalently for R=Rp. For the case α=3 the calculations are  very similar.(a) Hopf surfaceβ=βc and curves {l1=0}  (solid) and {l2=0} (dashed) for α=2.  (b) Hopf surface β=βc and curves {l1=0}  (solid)  and {l2=0}  (dashed) for α=3.
(a)(b)Figure 7
Bifurcation diagram of system (1.1) for α=2 at the point Q1.(a) Typical pointsV1, V2, V3, V4, V5, V6, V7, V8, and V9 in the parameter space near the point R1 (see Figure 7). (b) Phase portraits of system (1.1) for the flow restricted to the center manifold and its continuation at the points Vi,  i=1,…,9 of Figure 8(a).
(a)(b)Theorem 4.6.
Considerα=2 in system (1.1). Then there is a unique point Q1=Qα=2=(μ,ν,β2c) where the surfaces {l1=0} and {l2=0} intersect transversally on the Hopf hypersurface H2. For the parameter values at the point Q1, system (1.1) has a transversal Hopf point of codimension 3 at E+ which is unstable since l3(Q1)>0. The bifurcation diagram of system (1.1) at the point Q1 is illustrated in Figure 7. In Figure 8 the bifurcation diagram of system (1.1) at a typical point R1 of Figure 7 is drawn.Computer Assisted Proof
The pointQ1 is the intersection of the surfaces {l1=0} and {l2=0} on the Hopf hypersurface β=β2c taking into account α=2. It is defined and obtained as the solution of the equations g(μ,ν,α)=0 (see Theorem 4.3) and h(μ,ν,α)=0 (see Theorem 4.4). The point Q1 has coordinates
(4.26)μ=0.01528952519719747820…,ν=4.03634524181150646828…,β2c=4.02335387131139691747….
The existence and uniqueness of Q1 with the previous coordinates has been established numerically with the software MATHEMATICA 5 [18] (see also [17]) and can be checked by the Newton-Kantorovich Theorem [19].
For the pointQ1 take five decimal round-off coordinates μ=0.01528, ν=4.03634,  and β2c=4.02335 for simplicity. For these values of the parameters one has
(4.27)p=(-0.10542+0.46290i,0.21084+0.10811i,0.22356-0.43599i),q=(3.78949+0.96718i,0.48359i,4.02335),h11=(-371.10203,0,-371.10203),h20=(82.25866+60.23037i,-1.95710+26.99100i,112.27829+8.14125i),h30=(1209.86465+3423.09545i-366.52849+1469.19120i4074.40090+1016.46677i),G21=-2878.05753i,h21=(13925.76364+2288.71956i-623.13603-1297.11421i13152.94476+5184.29327i),h40=(-83816.60065+166465.53895i-47712.54683+79760.69093i165895.88420+99238.29710i),h31=(2.26342·106+1.12088·106i-87119.52533+376583.80975i2.56884·106+435080.93606i),h22=(-8.42248·106,0,-8.42248·106),G32=-9.92235·107i,h32=(5.84180·108+1.09060·108i-1.70175·107-6.11162·107i5.51875·108+2.34141·108i),h41=(1.25475·108+1.68063·108i-1.64672·107+6.39292·107i2.25787·108+5.77662·107i),h42=(1.93748·1011+9.38389·1010i-7.94067·109+2.60031·1010i2.15396·1011+4.67300·1010i),h33=(-7.28682·1011,0,-7.28682·1011),G43=1.17553·1010-1.14289·1013i.
From (3.9), (3.10), (3.11), and (4.27)  one has
(4.28)l1(Q1)=0,l2(Q1)=0,l3(Q1)=1144ReG43=1144(1.17553·1010).
The previous calculations have also been corroborated with 20 decimals round-off precision performed using the software MATHEMATICA 5 [18]; see [17].
The gradients of the functionsl1, given in (4.20), and l2, given in (4.23), at the point Q1 are, respectively,
(4.29)(181.95647,-14.19675),(-5.26376·106,4.12653·106).
The transversality condition at Q1 is equivalent to the nonvanishing of the determinant of the matrix whose columns are the above gradient vectors, which was evaluated and has furnished the value 6.76121·108>0. The transversality condition being satisfied, the bifurcation diagrams in Figures 7 and 8  follow from the work of Takens [16], taking into account the orientation and signs.
Definel0=β-β2c. The open region in the parameter space where system (1.1) has three small periodic orbits bifurcating from the equilibria E± can be described by  l2<0,  l1>0,  and l0>0 with |l2|≫l1≫l0>0. See Figure 7. The phase portraits of system (1.1) for the flow restricted to the center manifold and its continuation are shown in Figure 8(b). For parameters at  V1 the equilibrium is unstable; for parameters at  V2 the equilibrium is an weak unstable focus (Hopf point with positive l1); for parameters at  V3 the equilibrium is stable and an unstable limit cycle appears from the Hopf bifurcation; for parameters at V4 the equilibrium is an weak stable focus (Hopf point with negative l1) and there is an unstable limit cycle; for parameters at V5 the equilibrium is unstable and a stable limit cycle appears from the Hopf bifurcation, so there are two limit cycles encircling the equilibrium; for parameters at V6 the equilibrium is unstable and two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles; for parameters at V7 the equilibrium is stable and two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles; for parameters at V8 the equilibrium is stable and two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles; for parameters at V9 the equilibrium is stable and there are three limit cycles encircling the equilibrium.


## 4.1. Hopf Bifurcations atE0

In this subsection we study the stability ofE0 under the conditions β=β1c and -1<αμ<0, that is, on the Hopf hypersurface H1 complementary to the range of validity of Proposition 2.2.The Jacobian matrix of the functionf given in (2.1) at E0 is given by(4.1)A=Df(E0)=(μνν01-α-10β1c0).
Then, using the notation of the previous section (see expression (3.3)) the multilinear symmetric functions corresponding to f can be written as(4.2)B(x,y)=D(x,y,z,u)=E(x,y,z,u,v)=(0,0,0),C(x,y,z)=(-6νx1y1z1,0,0).The eigenvalues ofA are(4.3)λ1=-α+μν<0,λ2,3=±iω0=±iν-μ(1+αμ)α,
and from (3.5) one has
(4.4)q=(-iω0νμν-iω0,iω0,β1c),p̅=1ρ0(iω0(μν-iω0),-iω0(μν-iω0)2,(μν-iω0)2),
where
(4.5)ρ0=β1c(μν-iω0)2+(ν+(μν-iω0)2)ω02.Using these calculations we prove the next two theorems, from which statement (a) of the introduction follows.Theorem 4.1.
Consider the four-parameter family of differential equations (1.1). The first Lyapunov coefficient associated to the equilibrium E0 is given by
(4.6)l1(μ,ν,α,β1c)=-3α(αμ+1)ν3(α2-2μνα-ν)2(α3-2μνα2-μν2).
If
(4.7)s(μ,ν,α)=-α2+2μνα+ν
is different from zero then system (1.1) has a transversal Hopf point at E0 for β=β1c and -1<αμ<0.
More precisely, if(μ,ν,α,β1c)∈S∪U (see Figure 2(a)), then system (1.1) has a Hopf point of codimension 1 at E0. If (μ,ν,α,β1c)∈S then E0 is asymptotically stable and for each β<β1c, but close to β1c, there exists a stable periodic orbit near the unstable equilibrium point E0. If (μ,ν,α,β1c)∈U then E0 is unstable and for each β>β1c, but close to β1c, there exists an unstable periodic orbit near the asymptotically stable equilibrium point E0.Proof.
As the functionB(x,y)=(0,0,0) then the expression G21 in (3.9) reduces to G21〈p,C(q,q,q̅)〉. Performing the calculations we get
(4.8)G21=-6ν4ω04(μ2ν2+ω02)(β1c(μν-iω0)2+ω02(μ2ν2-2iμω0ν+ν-ω02)).
Substituting the expressions of ω0 and β1c into the above expression of G21 and taking the real part we have, from (3.9), the expression of the first Lyapunov coefficient l1 given in  (4.6). Note that the sign of l1 is given by the sign of the function s defined in (4.7).
It remains only to verify the transversality condition of the Hopf bifurcation. In order to do so, consider the family of differential equations (1.1) regarded as dependent on the parameter β. The real part, η=η(β), of the pair of complex eigenvalues at the critical parameter β=β1c verifies
(4.9)η′(β1c)=Re〈p,dAdβ|β=β1cq〉=-α22(α3-2μνα2-μν2)<0.
Therefore, the transversality condition at the Hopf point holds.
The setsS and U illustrated in Figure 2(a) are defined implicitly as {l1<0} and {l1>0},  respectively. The theorem is proved.Theorem 4.2.
Consider the four-parameter family of differential equations (1.1) restricted to β=β1c  and -1<αμ<0. Then the second Lyapunov coefficient l2 associated to the equilibrium E0 is negative on the Hopf hypersurface where the first Lyapunov coefficient l1 vanishes.Proof.
Following the notation introduced in Section3, we have B(x,y)=D(x,y,z,u)=E(x,y,z,u,v)=(0,0,0). Then
(4.10)ℋ32=3C(q,q,h̅21)+6C(q,q̅,h21)+3C(q̅,q̅,h30).
The expressions of the complex vectors h21 and h30 are too long to be put in print. The interested author can encounter such expressions and all the details of the calculations presented in this paper in the website [17].By direct calculations the second Lyapunov coefficientl2 associated to E0 can be written as(4.11)l2(μ,ν,α,β1c)=-3α2(αμ+1)2ν5N(μ,ν,α,β1c)2μ(-α3+2μνα2+μν2)3(-α3+2μνα2+8μ2ν2α+9μν2),
where N(μ,ν,α,β1c) has the form(4.12)α10-8μνα9+ν(16μ2ν-1)α8+2μν2(8μ2ν-7)α7+2μ2ν3(61-40μ2ν)α6+2μν3(32ν2μ4-112νμ2+13)α5-μ2ν4(64νμ2+19)α4+40μ3ν5(10μ2ν-9)α3+μ2ν5(-192ν2μ4+724νμ2-121)α2-8μ3ν6(43μ2ν-34)α-142μ4ν7.
Along the surface {l1=0} we have, from   (4.7) of Theorem 4.1, s(μ,ν,α)=-α2+2μνα+ν=0, or equivalently, να2/(1+2αμ). Substituting this expression into (4.11) one has(4.13)l2|{l1=0}=-9α11(αμ+1)(4αμ-1)(2αμ+1)4(8αμ-1)<0.Here the transversality condition is equivalent to the transversality of the hypersurfacesβ=β1c and  l1(μ,ν,α,β)=0. The bifurcation diagram for a typical point R on the surface {l1=0} is depicted in Figure 3.


## 4.2. Hopf Bifurcations atE±

In this subsection we study the stability ofE± under the conditions β=β2c and 0<2αμ<1, that is, on the Hopf hypersurface H2 complementary to the range of validity of Proposition 2.3.The Jacobian matrix of the functionf given in (2.1) calculated at E+ is given by(4.14)A=Df(E+)=(-2μνν01-α-10β2c0).
Then, using the notation of Section 3 (see expression (3.3)) the multilinear symmetric functions corresponding to f can be written as(4.15)B(x,y)=(-6νμx1y1,0,0),C(x,y,z)=(-6νx1y1z1,0,0),D=E=K=L=(0,0,0).The eigenvalues ofA are(4.16)λ1=-α-2μν<0,λ2,3=±iω0=±iν2μ(1-2αμ)α,
and from (3.5) one has(4.17)q=(iω0ν2μν+iω0,iω0,β2c),(4.18)p̅=1ρ0(-iω0(2μν+iω0),-iω0(2μν+iω0)2,(2μν+iω0)2),
where(4.19)ρ0=β2c(2μν+iω0)2+(ν+(2μν+iω0)2)ω02.Using these calculations we prove the next three theorems, from which statement (b) of introduction follows.Theorem 4.3.
Consider the four-parameter family of differential equations (1.1). The first Lyapunov coefficient associated to the equilibrium E+ is given by
(4.20)l1(μ,ν,α,β2c)=-3α(2αμ-1)ν3g(μ,ν,α)(α3+4μνα2+2μν2)(α3+4μνα2-12μ2ν2α+8μν2),
where
(4.21)g(μ,ν,α)=α5+8μνα4+ν(4μ2ν-1)α3+2μν2(5-24μ2ν)α2+80μ2ν3α-20μν3.
If g(μ,ν,α) is different from zero then system (1.1) has a transversal Hopf point at E+ for β=β2c and 0<2αμ<1.
More precisely, if(μ,ν,α,β2c)∈S∪U (see Figure 2(b)) then system (1.1) has a Hopf point of codimension 1 at E+. That is, if (μ,ν,α,β2c)∈S then E+ is asymptotically stable and for each β<β2c, but close to β2c, there exists a stable periodic orbit near the unstable equilibrium point E+; if (μ,ν,α,β2c)∈U then E+ is unstable and for each β>β2c, but close to β2c, there exists an unstable periodic orbit near the asymptotically stable equilibrium point E+.Proof.
The theorem follows by expanding the expressions in definition of the first Lyapunov coefficient (3.9). It relies on extensive calculations involving the vector q in (4.17), the vector p in (4.18) and the multilinear functions B and C. The inclusion here of the enormous expressions obtained from the mentioned calculations would not bring any contribution in clarifying the reading of the text. Then the detailed calculations related to this proof, corroborated by Computer Algebra, have been posted in [17], including its implementation using the software MATHEMATICA 5 [18].
Notice that the sign of the first Lyapunov coefficient is given by the sign of the function-g defined in (4.21) since the denominator of l1 and the expression 3(2αμ-1)ν3 are negative. In Figure 2(b) the surface {l1=0} and the regions S and U are illustrated. These regions are defined implicitly as {l1<0} and {l1>0}, respectively.
It remains only to verify the transversality condition of the Hopf bifurcation. In order to do so, consider the family of differential equations (1.1) regarded as dependent on the parameter β. The real part, η=η(β), of the pair of complex eigenvalues at the critical parameter β=β2c verifies
(4.22)η′(β2c)=Re〈p,dAdβ|β=β2cq〉=-α22(α3+4μνα2+2μν2)<0.
Therefore, the transversality condition holds at the Hopf point. The theorem is proved.Theorem 4.4.
Consider the four-parameter family of differential equations (1.1) restricted to β=β2c and 0<2αμ<1. Then the second Lyapunov coefficient l2 associated to the equilibrium E+ is given by
(4.23)l2(μ,ν,α,β2c)=α2(2αμ-1)ν5h(μ,ν,α)4μ(α3+4μνα2-32μ2ν2α+18μν2)d(μ,ν,α),
where the function h(μ,ν,α) is given by
(4.24)21μα20+(588μ2ν-4)α19+3μν(1876μ2ν-53)α18+4ν(3192ν2μ4-294νμ2+1)α17+μν2(-121296ν2μ4+14156νμ2-65)α16-12μ2ν3(59248ν2μ4-18704νμ2+473)α15+μν3(114240ν3μ6+716976ν2μ4-62932νμ2+199)α14+16μ2ν4(477120ν3μ6-165272ν2μ4+769νμ2+119)α13+4μ3ν5(1634304ν3μ6-4121136ν2μ4+751264νμ2-16945)α12-16μ2ν5(2225664ν4μ8-327552ν3μ6-585624ν2μ4+53587νμ2-262)α11-4μ3ν6(6967296ν4μ8-29565184ν3μ6+7648944ν2μ4+81592νμ2-11129)α10+16μ4ν7(4644864ν4μ8+773376ν3μ6-8268704ν2μ4+1627004νμ2-41405)α9+4μ3ν7(-80474112ν4μ8+32219456ν3μ6+13655344ν2μ4-1903712νμ2+17531)α8+16μ4ν8(5308416ν4μ8+34114176ν3μ6-13584144ν2μ4-51152νμ2+39713)α7+32μ5ν9(1327104ν4μ8-10505088ν3μ6-13899000ν2μ4+4382138νμ2-130059)α6-16μ4ν9(11446272ν4μ8-35076672ν3μ6-10290192ν2μ4+2504248νμ2-37769)α5+32μ5ν10(10175616ν3μ6-15588272ν2μ4-317442νμ2+129593)α4-64μ6ν11(4720704ν2μ4-3846404νμ2+142349)α3+128μ5ν11(1199432ν2μ4-495365νμ2+13321)α2-256μ6ν12(157944μ2ν-25841)α+4308992μ7ν13,
and the function d(μ,ν,α) in the denominator is given by
(4.25)(α6+8μνα5+4μ2ν2α4+2μν2(5-24μ2ν)α3+40μ2ν3α2-24μ3ν4α+16μ2ν4)3.
If (μ,ν,α,β2c)∈C1∪C2 (see Figure 2(b)), then system (1.1) has a transversal Hopf point of codimension 2 at E+. More specifically, if (μ,ν,α,β2c)∈C1 then the Hopf point at E+ is asymptotically stable and the bifurcation diagram is illustrated in Figure 3 for a typical point R. If (μ,ν,α,β2c)∈C2 then the Hopf point at E+ is unstable and the bifurcation diagram is drawn in Figure 4 for a typical point T.Proof.
The theorem follows by expanding the expressions in definition (3.10) of the second Lyapunov coefficient. It relies on extensive calculation involving the vector q in (4.17), the vector p in (4.18) and the multilinear functions B, C, D, and E. The calculations in this proof, corroborated by Computer Algebra, have been posted in [17].
In Figure2 we present a geometric synthesis interpreting the long calculations involved in this proof. The sign of -h(μ,ν,α) gives the sign of the second Lyapunov coefficient since the function d(α,μ,ν) in the denominator of l2 is positive. The graph of h(μ,ν,α)=0 on the surface {l1=0}, which determines the curve C, and the signs of the first and second Lyapunov coefficients are also illustrated. As follows, l2<0 on the open region on the surface {l1=0}, denoted by C1. On this region a typical reference point R is depicted. Also l2>0 on the open region on the surface {l1=0}, denoted by C2. This region contains the typical reference point, denoted by T.
The bifurcation diagrams of system (1.1) at the points R and T are illustrated in Figures 3 and 4, respectively. The theorem is proved.Remark 4.5.
There are numerical evidences that the sign of the third Lyapunov coefficient is always positive along the curveC given by the intersection of the surfaces {l1=0} and {l2=0} (see Figure 2(b)). For α=1 see the article in [5]. For α=2 see Theorem 4.6 in what follows. Figure 5 shows the behavior of the third Lyapunov coefficient as a function of α, for (μ,ν,β) on the curve C. So the bifurcation diagrams depicted in Figures 7 and 8 are essentially the same for all values of α≥1.Figure 5
Graph ofl3 with respect to α.All the calculations presented here are illustrated in Figure6 for the cases α=2 (a) and α=3 (b). In this figure are depicted the surfaces β=β2c, the solid curves where l1=0, the regions S and U where l1<0 and l1>0, respectively, the dashed curves where l2=0 and the arcs C1 and C2 where l2<0 and l2>0, respectively. The bifurcation diagrams for typical points R1 and R2 are depicted in Figure 3 while the bifurcation diagrams for typical points T1 and T2 are shown in Figure 4. In the next theorem we analyze the sign of the third Lyapunov coefficient at the point Q1, that is, for the case α=2 or equivalently for R=Rp. For the case α=3 the calculations are  very similar.(a) Hopf surfaceβ=βc and curves {l1=0}  (solid) and {l2=0} (dashed) for α=2.  (b) Hopf surface β=βc and curves {l1=0}  (solid)  and {l2=0}  (dashed) for α=3.
(a)(b)Figure 7
Bifurcation diagram of system (1.1) for α=2 at the point Q1.(a) Typical pointsV1, V2, V3, V4, V5, V6, V7, V8, and V9 in the parameter space near the point R1 (see Figure 7). (b) Phase portraits of system (1.1) for the flow restricted to the center manifold and its continuation at the points Vi,  i=1,…,9 of Figure 8(a).
(a)(b)Theorem 4.6.
Considerα=2 in system (1.1). Then there is a unique point Q1=Qα=2=(μ,ν,β2c) where the surfaces {l1=0} and {l2=0} intersect transversally on the Hopf hypersurface H2. For the parameter values at the point Q1, system (1.1) has a transversal Hopf point of codimension 3 at E+ which is unstable since l3(Q1)>0. The bifurcation diagram of system (1.1) at the point Q1 is illustrated in Figure 7. In Figure 8 the bifurcation diagram of system (1.1) at a typical point R1 of Figure 7 is drawn.Computer Assisted Proof
The pointQ1 is the intersection of the surfaces {l1=0} and {l2=0} on the Hopf hypersurface β=β2c taking into account α=2. It is defined and obtained as the solution of the equations g(μ,ν,α)=0 (see Theorem 4.3) and h(μ,ν,α)=0 (see Theorem 4.4). The point Q1 has coordinates
(4.26)μ=0.01528952519719747820…,ν=4.03634524181150646828…,β2c=4.02335387131139691747….
The existence and uniqueness of Q1 with the previous coordinates has been established numerically with the software MATHEMATICA 5 [18] (see also [17]) and can be checked by the Newton-Kantorovich Theorem [19].
For the pointQ1 take five decimal round-off coordinates μ=0.01528, ν=4.03634,  and β2c=4.02335 for simplicity. For these values of the parameters one has
(4.27)p=(-0.10542+0.46290i,0.21084+0.10811i,0.22356-0.43599i),q=(3.78949+0.96718i,0.48359i,4.02335),h11=(-371.10203,0,-371.10203),h20=(82.25866+60.23037i,-1.95710+26.99100i,112.27829+8.14125i),h30=(1209.86465+3423.09545i-366.52849+1469.19120i4074.40090+1016.46677i),G21=-2878.05753i,h21=(13925.76364+2288.71956i-623.13603-1297.11421i13152.94476+5184.29327i),h40=(-83816.60065+166465.53895i-47712.54683+79760.69093i165895.88420+99238.29710i),h31=(2.26342·106+1.12088·106i-87119.52533+376583.80975i2.56884·106+435080.93606i),h22=(-8.42248·106,0,-8.42248·106),G32=-9.92235·107i,h32=(5.84180·108+1.09060·108i-1.70175·107-6.11162·107i5.51875·108+2.34141·108i),h41=(1.25475·108+1.68063·108i-1.64672·107+6.39292·107i2.25787·108+5.77662·107i),h42=(1.93748·1011+9.38389·1010i-7.94067·109+2.60031·1010i2.15396·1011+4.67300·1010i),h33=(-7.28682·1011,0,-7.28682·1011),G43=1.17553·1010-1.14289·1013i.
From (3.9), (3.10), (3.11), and (4.27)  one has
(4.28)l1(Q1)=0,l2(Q1)=0,l3(Q1)=1144ReG43=1144(1.17553·1010).
The previous calculations have also been corroborated with 20 decimals round-off precision performed using the software MATHEMATICA 5 [18]; see [17].
The gradients of the functionsl1, given in (4.20), and l2, given in (4.23), at the point Q1 are, respectively,
(4.29)(181.95647,-14.19675),(-5.26376·106,4.12653·106).
The transversality condition at Q1 is equivalent to the nonvanishing of the determinant of the matrix whose columns are the above gradient vectors, which was evaluated and has furnished the value 6.76121·108>0. The transversality condition being satisfied, the bifurcation diagrams in Figures 7 and 8  follow from the work of Takens [16], taking into account the orientation and signs.
Definel0=β-β2c. The open region in the parameter space where system (1.1) has three small periodic orbits bifurcating from the equilibria E± can be described by  l2<0,  l1>0,  and l0>0 with |l2|≫l1≫l0>0. See Figure 7. The phase portraits of system (1.1) for the flow restricted to the center manifold and its continuation are shown in Figure 8(b). For parameters at  V1 the equilibrium is unstable; for parameters at  V2 the equilibrium is an weak unstable focus (Hopf point with positive l1); for parameters at  V3 the equilibrium is stable and an unstable limit cycle appears from the Hopf bifurcation; for parameters at V4 the equilibrium is an weak stable focus (Hopf point with negative l1) and there is an unstable limit cycle; for parameters at V5 the equilibrium is unstable and a stable limit cycle appears from the Hopf bifurcation, so there are two limit cycles encircling the equilibrium; for parameters at V6 the equilibrium is unstable and two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles; for parameters at V7 the equilibrium is stable and two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles; for parameters at V8 the equilibrium is stable and two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles; for parameters at V9 the equilibrium is stable and there are three limit cycles encircling the equilibrium.


## 5. Numerical Simulations

In this section we present some numerical simulations of system (1.1) for several values of the parameters ν, μ,  α, and β. The main purpose is to illustrate the creation of stable and unstable limit cycles through the Hopf bifurcations at the equilibria E0 and E±, proved to occur in the previous sections. The simulations were developed using the Software MAPLE 8, under the Runge-Kutta method of fourth order with several different stepsizes. We also perform a numerical study of the continuation of periodic orbits which arise from the Hopf bifurcations at the point E±. Throughout this analysis we have found period doubling and homoclinic bifurcations which seems to lead to the creation of strange attractors for (1.1), for some parameter values. For the study of homoclinic bifurcations in three-dimensional systems see [20, 21].


### 5.1. Bifurcations atE0

Forμ<0, system (1.1) has the  origin E0=(0,0,0) as its unique  equilibrium point. According to Theorem 4.1, the system undergoes a Hopf bifurcation when the parameter β crosses the critical value β1c=ν(α-μν)(1+αμ)/α with -1<αμ<0. This type of bifurcation is illustrated in Figures 9 and 10. To draw these figures we have taken ν=3, μ=-0.25, α=2 and have varied the parameter β, whose values are presented in the captions of the figures. Observe that for this parameter values we have β1c=2.0625. One can observe that when the parameter β tends to the critical value β=β1c from above the spiralling behavior of the solution becomes slower as expected, corroborating the correctness of the calculations presented in the previous sections (see Figures 9(a) and 9(b)); the limit cycle which arises when the parameter β crosses β=β1c is very small (observe the range of the state variables (x,y,z) in Figure 10). The numerical analysis performed suggests that such small limit cycle which birth in the Hopf bifurcation for β=β1c  dies for β=0; that is, the cycle exists for β∈(0,2.0625). One can observe that as the parameter β stay away from the critical value β1c the solutions tend faster to the limit cycle (see Figure 11). The shape and behavior of the limit cycles in Figure 11 resemble the one of the Van der Pol system, see [22, page 267].Phase portrait of system (1.1) near the equilibrium point E0 for the parameter values ν=3, μ=-0.25, α=2, and  (a) β=β1c+0.5=2.5625; (b) β=β1c+0.05=2.1125. Time of integration is [5,350] for (a) [5,450] for (b). Initial condition is (0.001,-0.001,0.001). Stepsize = 0.1. In both cases the point E0 is a stable focus, but the convergence in the case (b) is slower than in (a) since the parameter β is closer to the critical value β1c=2.0625. Eigenvalues at E0 are (a) λ1=-2.624, λ2,3=-0.062±0.853i; (b) λ1=-2.737, λ2,3=-0.006±0.760i.
(a)(b)Limit cycle created in the Hopf bifurcation for system (1.1) at the point E0. The parameter values are  ν=3, μ=-0.25, α=2, and β=β1c-0.1=1.9625. Time of integration is  [450,750] in (a) and [900,950] in (b). Initial condition is (0.001,-0.001,0.001). Stepsize = 0.1. The point E0 is an unstable focus and a limit cycle arises around it. Eigenvalues at E0 are λ1=-2.774,λ2,3=0.012±0.728i.
(a)(b)Limit cycle created in the Hopf bifurcation for system (1.1) at the point E0. The parameter values are ν=3, μ=-0.25, α=2, and  β=β1c-1=1.0625 in (a); β=β1c-2=0.0625 in (b). Time of integration is [0,150]  in  (a) and [0,100] in (b). Initial condition is (0.001,-0.001,0.001). Stepsize = 0.1. As the parameter β stay away from the critical value β1c the solutions tend faster to the limit cycle. Eigenvalues at E0  are  (a) λ1=-2.985, λ2,3=0.117±0.502i; (b) λ1=-3.203, λ2=0.418, λ3=0.034.
(a)(b)


### 5.2. Bifurcations atE±

Forμ>0, system (1.1) has the origin E0 and also the symmetric equilibria E±=(±μ,0,±μ). According to Theorem 4.3, the system undergoes a Hopf bifurcation when the parameter β crosses the critical value β=β2c=ν(α+2μν)(1-2αμ)/α, for 0<αμ<1/2. This type of bifurcation is illustrated in Figures 12 and 13, where we have taken ν=10, μ=0.1, α=2 and have varied the parameter β, whose values are presented in the captions of the figures. Observe that for these values of the parameters ν, μ, and α we have β2c=12. In the figures are shown the 1-dimensional unstable manifolds of E0, which spiral to the equilibria E± for β>β2c. Observe that when the parameter β tends to the critical value β=β2c the spiralling behavior of the solutions  becomes slower (Figure 12(b)); for β<β2c the unstable manifolds spiral to the symmetric limit cycles   arise in the Hopf bifurcation for β=β2c (Figure 13).Phase portraits of system (1.1) before the Hopf bifurcation at the points E±. The parameter values are ν=10, μ=0.1, α=2, and   β=β2c+0.8=12.8  in (a); β=β2c+0.2=12.2  in (b). Time of integration is [-4.5,50] for (a); [-4.5,270] for (b). Initial conditions are (0.0001,0,0.0001) and (-0.0001,0,-0.0001). Stepsize = 0.05. Eigenvalues at E0 are (a) λ1=1.951, λ2,3=-1.475±2.093i; (b) λ1=1.987, λ2,3=-1.493±1.976i. Eigenvalues at E± are (a) λ1=-3.928, λ2,3=-0.035±2.552i; (b) λ1=-3.981, λ2,3=-0.009±2.475i.
(a)(b)Phase portraits of system (1.1) after the Hopf  bifurcation at the points E±. The parameter values are  ν=10, μ=0.1, α=2, and  β=β2c-0.2=11.8. The equilibria become unstable and two symmetric limit cycles arise. Time of integration is [-4.5,270] for (a);  [500,550] for (b). Initial conditions are (0.0001,0,0.0001),  (-0.0001,0,-0.0001),  (0.5,0.01,0.5), and (-0.5,0.01,-0.5) for (a); (0.0001,0,0.0001) and (-0.0001,0,-0.0001) for (b). Stepsize = 0.05. Eigenvalues at E0 are  λ1=2.012, λ2,3=-1.506±1.895i. Eigenvalues at E±  are  λ1=-4.018, λ2,3=0.009±2.423i.
(a)(b)We observe that  when the parameterβ moves away from the critical value β2c, the numerical simulations performed suggest that  firstly a duplication of period occurs with the limit cycles and then two double symmetric homoclinic loops to the origin E0 appear (see Figure 14). After this, a strange attractor seems to exist near these homoclinic  loops (Figure 15). This is one of the mechanisms through which system (1.1) enters into chaotic regime. Observe that it begins with the creation of the limit cycles in the Hopf bifurcations which takes place at the points E± for the critical value of parameter β=β2c. For the case α=1 another mechanism giving rise to chaotic attractor is described in [5].(a) Period doubling bifurcation for the limit cycle created in the Hopf bifurcation at the pointE±. The parameter values are  ν=10, μ=0.1, α=2, and  β=β2c-1.13=10.87. (b) Double symmetric homoclinic loop to the origin. Parameter values ν=10, μ=0.1, α=2, and  β=β2c-1.214322=10.785678. Initial conditions are (0.0001,0,0.0001) and (-0.0001,0,-0.0001). Eigenvalues at E0 are (a) λ1=2.073, λ2,3=-1.536±1.697i; (b) λ1=2.079, λ2,3=-1.539±1.678i. Eigenvalues at E± are  (a) λ1=-4.104, λ2,3=0.052±2.300i; (b) λ1=-4.111, λ2,3=0.055±2.289i.
(a)(b)Figure 15
Strange attractor created after the existence of the double homoclinic loop shown in Figure14. Parameter values  are ν=10, μ=0.1, α=2, and  β=β2c-1.743=10.257. Initial conditions  are (0.0001,0,0.0001) and (-0.0001,0,-0.0001). Eigenvalues at E0 are  λ1=2.115, λ2,3=-1.557±1.556i; eigenvalues at E± are λ1=-4.161, λ2,3=0.080±2.218i.


### 5.3. Bifurcation at the Critical PointQ1 in the Parameter Space

As stated in Theorem4.6, the most complex Hopf bifurcation at the points E± occurs at the point Q1=(μ,ν,β2c) in the parameter space, when we consider α=2 fixed. In this subsection we present some numerical simulations of the solutions of system (1.1) for (μ,ν,β) near and at the point Q1=(μQ1,νQ1,βQ1)=(0.015289525,4.036345241,4.023353871). In Figure 16 we show an approximation of the unstable manifolds of the equilibrium E0. These manifolds spiral to the equilibria E± if β>βQ1=4.023353871 and this spiralling behavior become slower as β tends to βQ1=β2c from above (for μ=μQ1 and ν=νQ1 fixed).Phase portraits of system (1.1) before the Hopf bifurcation at the points E±. The parameter values are near the point Q1 of Theorem 4.6 in the parameter space: μ=μQ1=0.015289525, ν=νQ1=4.036345241,  and   β=βQ+0.5=4.523353871 in (a);  β=βQ+0.01=4.033353871 in (b). Time of integration is [-4,50] (a);  [-4,1380]  (b). Initial conditions are (0.001,0,0.001) and (-0.001,0,-0.001). Stepsize = 0.01. Eigenvalues at E0 are (a) λ1=0.280, λ2=-1.602, λ3=-0.617; (b) λ1=0.357, λ2=-1.937, λ3=-0.358. Eigenvalues at E± are (a) λ1=-1.891, λ2,3=-0.115±0.530i; (b) λ1=-2.119, λ2,3=-0.002±0.484i.
(a)(b)For(μ,α,β) at the point Q1 (see Theorem 4.6) the system presents a highly degenerate behavior. In fact, due to the vanishing of the first and second Lyapunov coefficients, l1 and l2, the solutions near the equilibria E± behave like periodic solutions, for a long time of integration. This behavior is shown in Figure 17. Although we have proved theoretically the instability of the equilibria E± in Theorem 4.6, since l3(Q1)>0, the solutions move away from these points very slowly, with repeller rate as 10-9. Then, it is not detected numerically  and the solutions near the equilibria E± seem to be periodic as shown in Figure 17. In addition there is also a “big" periodic orbit encircling the equilibria E±. It would be very interesting to observe how the real electronic circuit behave for the parameter values at the bifurcation point Q1.Figure 17
Phase portrait of system (1.1) for the parameters at the point Q1=(μQ1,νQ1,βQ1) of Theorem 4.6 in the parameter space: (μ,ν.β). The solutions near the equilibria E± behave like periodic solutions. Time of integration is [10500,12500]. There are several initial conditions. Stepsize = 0.1. There is also a “big” periodic orbit encircling the equilibria E±. Eigenvalues at E0 are  λ1=0.359,λ2=-1.942,λ3=-0.355. Eigenvalues at E±  are λ1=-2.123, λ2,3=-0.312×10-10±0.483i≈±0.483i.Forβ<βQ1, the equilibria become unstable and the solutions starting near them tend to a “big” periodic orbit (see Figure 18), which seems to be the unique attractor of the system in this case.Phase portrait of system (1.1) for μ=μQ1, ν=νQ1, and β=3.973353871<βQ1. The three equilibria are unstable and there exists a “big” stable periodic orbit encircling them. Eigenvalues at E0  are  λ1=0.368, λ2=-1.969, λ3=-0.337. Eigenvalues at E± are  λ1=-2.144, λ2,3=0.010±0.478i.
(a)(b)


## 5.1. Bifurcations atE0

Forμ<0, system (1.1) has the  origin E0=(0,0,0) as its unique  equilibrium point. According to Theorem 4.1, the system undergoes a Hopf bifurcation when the parameter β crosses the critical value β1c=ν(α-μν)(1+αμ)/α with -1<αμ<0. This type of bifurcation is illustrated in Figures 9 and 10. To draw these figures we have taken ν=3, μ=-0.25, α=2 and have varied the parameter β, whose values are presented in the captions of the figures. Observe that for this parameter values we have β1c=2.0625. One can observe that when the parameter β tends to the critical value β=β1c from above the spiralling behavior of the solution becomes slower as expected, corroborating the correctness of the calculations presented in the previous sections (see Figures 9(a) and 9(b)); the limit cycle which arises when the parameter β crosses β=β1c is very small (observe the range of the state variables (x,y,z) in Figure 10). The numerical analysis performed suggests that such small limit cycle which birth in the Hopf bifurcation for β=β1c  dies for β=0; that is, the cycle exists for β∈(0,2.0625). One can observe that as the parameter β stay away from the critical value β1c the solutions tend faster to the limit cycle (see Figure 11). The shape and behavior of the limit cycles in Figure 11 resemble the one of the Van der Pol system, see [22, page 267].Phase portrait of system (1.1) near the equilibrium point E0 for the parameter values ν=3, μ=-0.25, α=2, and  (a) β=β1c+0.5=2.5625; (b) β=β1c+0.05=2.1125. Time of integration is [5,350] for (a) [5,450] for (b). Initial condition is (0.001,-0.001,0.001). Stepsize = 0.1. In both cases the point E0 is a stable focus, but the convergence in the case (b) is slower than in (a) since the parameter β is closer to the critical value β1c=2.0625. Eigenvalues at E0 are (a) λ1=-2.624, λ2,3=-0.062±0.853i; (b) λ1=-2.737, λ2,3=-0.006±0.760i.
(a)(b)Limit cycle created in the Hopf bifurcation for system (1.1) at the point E0. The parameter values are  ν=3, μ=-0.25, α=2, and β=β1c-0.1=1.9625. Time of integration is  [450,750] in (a) and [900,950] in (b). Initial condition is (0.001,-0.001,0.001). Stepsize = 0.1. The point E0 is an unstable focus and a limit cycle arises around it. Eigenvalues at E0 are λ1=-2.774,λ2,3=0.012±0.728i.
(a)(b)Limit cycle created in the Hopf bifurcation for system (1.1) at the point E0. The parameter values are ν=3, μ=-0.25, α=2, and  β=β1c-1=1.0625 in (a); β=β1c-2=0.0625 in (b). Time of integration is [0,150]  in  (a) and [0,100] in (b). Initial condition is (0.001,-0.001,0.001). Stepsize = 0.1. As the parameter β stay away from the critical value β1c the solutions tend faster to the limit cycle. Eigenvalues at E0  are  (a) λ1=-2.985, λ2,3=0.117±0.502i; (b) λ1=-3.203, λ2=0.418, λ3=0.034.
(a)(b)


## 5.2. Bifurcations atE±

Forμ>0, system (1.1) has the origin E0 and also the symmetric equilibria E±=(±μ,0,±μ). According to Theorem 4.3, the system undergoes a Hopf bifurcation when the parameter β crosses the critical value β=β2c=ν(α+2μν)(1-2αμ)/α, for 0<αμ<1/2. This type of bifurcation is illustrated in Figures 12 and 13, where we have taken ν=10, μ=0.1, α=2 and have varied the parameter β, whose values are presented in the captions of the figures. Observe that for these values of the parameters ν, μ, and α we have β2c=12. In the figures are shown the 1-dimensional unstable manifolds of E0, which spiral to the equilibria E± for β>β2c. Observe that when the parameter β tends to the critical value β=β2c the spiralling behavior of the solutions  becomes slower (Figure 12(b)); for β<β2c the unstable manifolds spiral to the symmetric limit cycles   arise in the Hopf bifurcation for β=β2c (Figure 13).Phase portraits of system (1.1) before the Hopf bifurcation at the points E±. The parameter values are ν=10, μ=0.1, α=2, and   β=β2c+0.8=12.8  in (a); β=β2c+0.2=12.2  in (b). Time of integration is [-4.5,50] for (a); [-4.5,270] for (b). Initial conditions are (0.0001,0,0.0001) and (-0.0001,0,-0.0001). Stepsize = 0.05. Eigenvalues at E0 are (a) λ1=1.951, λ2,3=-1.475±2.093i; (b) λ1=1.987, λ2,3=-1.493±1.976i. Eigenvalues at E± are (a) λ1=-3.928, λ2,3=-0.035±2.552i; (b) λ1=-3.981, λ2,3=-0.009±2.475i.
(a)(b)Phase portraits of system (1.1) after the Hopf  bifurcation at the points E±. The parameter values are  ν=10, μ=0.1, α=2, and  β=β2c-0.2=11.8. The equilibria become unstable and two symmetric limit cycles arise. Time of integration is [-4.5,270] for (a);  [500,550] for (b). Initial conditions are (0.0001,0,0.0001),  (-0.0001,0,-0.0001),  (0.5,0.01,0.5), and (-0.5,0.01,-0.5) for (a); (0.0001,0,0.0001) and (-0.0001,0,-0.0001) for (b). Stepsize = 0.05. Eigenvalues at E0 are  λ1=2.012, λ2,3=-1.506±1.895i. Eigenvalues at E±  are  λ1=-4.018, λ2,3=0.009±2.423i.
(a)(b)We observe that  when the parameterβ moves away from the critical value β2c, the numerical simulations performed suggest that  firstly a duplication of period occurs with the limit cycles and then two double symmetric homoclinic loops to the origin E0 appear (see Figure 14). After this, a strange attractor seems to exist near these homoclinic  loops (Figure 15). This is one of the mechanisms through which system (1.1) enters into chaotic regime. Observe that it begins with the creation of the limit cycles in the Hopf bifurcations which takes place at the points E± for the critical value of parameter β=β2c. For the case α=1 another mechanism giving rise to chaotic attractor is described in [5].(a) Period doubling bifurcation for the limit cycle created in the Hopf bifurcation at the pointE±. The parameter values are  ν=10, μ=0.1, α=2, and  β=β2c-1.13=10.87. (b) Double symmetric homoclinic loop to the origin. Parameter values ν=10, μ=0.1, α=2, and  β=β2c-1.214322=10.785678. Initial conditions are (0.0001,0,0.0001) and (-0.0001,0,-0.0001). Eigenvalues at E0 are (a) λ1=2.073, λ2,3=-1.536±1.697i; (b) λ1=2.079, λ2,3=-1.539±1.678i. Eigenvalues at E± are  (a) λ1=-4.104, λ2,3=0.052±2.300i; (b) λ1=-4.111, λ2,3=0.055±2.289i.
(a)(b)Figure 15
Strange attractor created after the existence of the double homoclinic loop shown in Figure14. Parameter values  are ν=10, μ=0.1, α=2, and  β=β2c-1.743=10.257. Initial conditions  are (0.0001,0,0.0001) and (-0.0001,0,-0.0001). Eigenvalues at E0 are  λ1=2.115, λ2,3=-1.557±1.556i; eigenvalues at E± are λ1=-4.161, λ2,3=0.080±2.218i.


## 5.3. Bifurcation at the Critical PointQ1 in the Parameter Space

As stated in Theorem4.6, the most complex Hopf bifurcation at the points E± occurs at the point Q1=(μ,ν,β2c) in the parameter space, when we consider α=2 fixed. In this subsection we present some numerical simulations of the solutions of system (1.1) for (μ,ν,β) near and at the point Q1=(μQ1,νQ1,βQ1)=(0.015289525,4.036345241,4.023353871). In Figure 16 we show an approximation of the unstable manifolds of the equilibrium E0. These manifolds spiral to the equilibria E± if β>βQ1=4.023353871 and this spiralling behavior become slower as β tends to βQ1=β2c from above (for μ=μQ1 and ν=νQ1 fixed).Phase portraits of system (1.1) before the Hopf bifurcation at the points E±. The parameter values are near the point Q1 of Theorem 4.6 in the parameter space: μ=μQ1=0.015289525, ν=νQ1=4.036345241,  and   β=βQ+0.5=4.523353871 in (a);  β=βQ+0.01=4.033353871 in (b). Time of integration is [-4,50] (a);  [-4,1380]  (b). Initial conditions are (0.001,0,0.001) and (-0.001,0,-0.001). Stepsize = 0.01. Eigenvalues at E0 are (a) λ1=0.280, λ2=-1.602, λ3=-0.617; (b) λ1=0.357, λ2=-1.937, λ3=-0.358. Eigenvalues at E± are (a) λ1=-1.891, λ2,3=-0.115±0.530i; (b) λ1=-2.119, λ2,3=-0.002±0.484i.
(a)(b)For(μ,α,β) at the point Q1 (see Theorem 4.6) the system presents a highly degenerate behavior. In fact, due to the vanishing of the first and second Lyapunov coefficients, l1 and l2, the solutions near the equilibria E± behave like periodic solutions, for a long time of integration. This behavior is shown in Figure 17. Although we have proved theoretically the instability of the equilibria E± in Theorem 4.6, since l3(Q1)>0, the solutions move away from these points very slowly, with repeller rate as 10-9. Then, it is not detected numerically  and the solutions near the equilibria E± seem to be periodic as shown in Figure 17. In addition there is also a “big" periodic orbit encircling the equilibria E±. It would be very interesting to observe how the real electronic circuit behave for the parameter values at the bifurcation point Q1.Figure 17
Phase portrait of system (1.1) for the parameters at the point Q1=(μQ1,νQ1,βQ1) of Theorem 4.6 in the parameter space: (μ,ν.β). The solutions near the equilibria E± behave like periodic solutions. Time of integration is [10500,12500]. There are several initial conditions. Stepsize = 0.1. There is also a “big” periodic orbit encircling the equilibria E±. Eigenvalues at E0 are  λ1=0.359,λ2=-1.942,λ3=-0.355. Eigenvalues at E±  are λ1=-2.123, λ2,3=-0.312×10-10±0.483i≈±0.483i.Forβ<βQ1, the equilibria become unstable and the solutions starting near them tend to a “big” periodic orbit (see Figure 18), which seems to be the unique attractor of the system in this case.Phase portrait of system (1.1) for μ=μQ1, ν=νQ1, and β=3.973353871<βQ1. The three equilibria are unstable and there exists a “big” stable periodic orbit encircling them. Eigenvalues at E0  are  λ1=0.368, λ2=-1.969, λ3=-0.337. Eigenvalues at E± are  λ1=-2.144, λ2,3=0.010±0.478i.
(a)(b)


## 6. Concluding Remarks

In this paper we analyze the Lyapunov stability of the equilibriaE0 and E± of a Van der Pol-Duffing system with a parallel resistor given by (1.1). Through these analyses we obtain the hypersurfaces in the 4-dimensional parameter space for which the system presents Hopf bifurcations at these equilibria (Figure 2). Then we make an extension of the analysis to the more degenerate cases, happening in the locus on the Hopf hypersurfaces where the Lyapunov coefficients associated to the Hopf bifurcations vanish. The bifurcation analysis at E0 (resp., E±) is pushed forward to the calculation of the second (resp., second and third) Lyapunov coefficients, which makes possible the determination of the Lyapunov stability at the equilibria as well as the largest number of small periodic orbits bifurcating of these points.With the analytic data provided in the analysis performed here, the bifurcation diagrams are established. Figures3, 4, and 8 provide a qualitative synthesis of the dynamical conclusions achieved for the parameter values where system (1.1) has the most complex equilibria. For parameters inside the open “tongue" (see Figures 7 and 8) there are one attracting periodic orbit and two other unstable periodic orbits coexisting with the attracting equilibrium points E±.The calculations of the Lyapunov coefficients, being extensive, rely on Computer Algebra and numerical evaluations carried out with the software MATHEMATICA 5 [18]. In the website [17] have been posted the main steps of the calculations in the form of notebooks for MATHEMATICA 5 as well as in pdf format.Moreover numerical simulations were performed for several values of the parameters, which illustrate and corroborate some of the analytical results stated. The most significant phase portraits are shown in Section6. Some global phenomena suggested by the numerical analysis performed, as the existence of strange attractors and a double symmetric homoclinic loop for system (1.1), are also briefly commented in Section 6. Although it is out of the main purpose of this note, we have included it here because the existence of these attractors is in some sense related to the Hopf bifurcations which occur at the equilibria.



---

*Source: 149563-2009-11-19.xml*