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

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.

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) where (𝑥,𝑦,𝑧)∈ℝ³ 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) by the following changes in variables, parameters, and rescaling in time:

𝑥=√𝑏𝑅𝑣₁,𝑦=√𝑏𝑅𝑣₂,𝑧=𝑏𝑅³𝑖_𝐿,𝜇=−(1+𝑎𝑅),𝜈=𝑐₂/𝑐₁,𝛽=𝑐₂𝑅²/𝐿,𝜏=𝑡/𝑅𝑐₂. (1.3)

In (1.2) 𝑣₁ and 𝑣₂ are the voltages across the capacitors 𝐶₁ and 𝐶₂, respectively, 𝑖_𝑁 is the current through the nonlinear diode 𝑁, 𝑖_𝐿 is the current through the inductor, 𝑐₁ and 𝑐₂ are the capacitances, 𝐿 is the inductance, 𝑅 is the linear resistor, 𝑅_𝑝 is the parallel resistor, 𝛼=1+(𝑅/𝑅_𝑝), and the prime denotes derivatives with respect to the time 𝑡. It is assumed that the current through the nonlinear diode is given by 𝑖_𝑁=𝑖_𝑁(𝑣₁)=𝑎𝑣₁+𝑏𝑣³₁, with 𝑎<0 and 𝑏>0. See Figure 1 and [1] for more details.

Figure 1 

ADVP circuit with parallel resistor 𝑅𝑝.

Despite the simplicity of the electronic circuit shown in Figure 1, 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 𝐸₀=(0,0,0), which exists for any parameter values. For 𝜇>0 it also has the symmetric equilibria 𝐸_±=(±√𝜇,0,±√𝜇). Codimension one Hopf bifurcation which occurs at the equilibrium point 𝐸₀ 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 equilibria 𝐸₀ and 𝐸_± 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 equilibrium 𝐸₀=(0,0,0) the Hopf hypersurface 𝐻₁ is obtained in the space of parameters (𝜇,𝜈,𝛼,𝛽) and the first Lyapunov coefficient 𝑙₁ is calculated. It is shown that this coefficient vanishes on a 2-dimensional surface contained in 𝐻₁, giving rise to codimension two Hopf bifurcations (see Figure 2(a)). Then the second Lyapunov coefficient 𝑙₂ is calculated and it is established that this coefficient is always negative on the surface {𝑙₁=0}.(b)For the symmetric equilibria 𝐸_±=(±√𝜇,0,±√𝜇) the Hopf hypersurface 𝐻₂ is obtained in the space of parameters (𝜇,𝜈,𝛼,𝛽), the first Lyapunov coefficient 𝑙₁ is calculated, and it is shown that this coefficient vanishes on a 2-dimensional surface contained in 𝐻₂ (see Figure 2(b)), giving rise to codimension two bifurcations. The second Lyapunov coefficient 𝑙₂ is calculated and it is established that this coefficient also vanishes on a 2-dimensional surface contained in 𝐻₂ (see Figure 2(b)), giving rise to codimension three bifurcations. The third Lyapunov coefficients along the curve 𝐶 given by the intersection of the surfaces {𝑙₁=0} and {𝑙₂=0} are obtained and found to be positive.

(a)

(b)

Figure 2 

Surfaces {𝑙1=0} for 𝐸0 (a) and 𝐸± (b), 𝐶={𝑙1=0}∩{𝑙2=0}.

The corresponding bifurcation diagrams are traced for the above bifurcations (see Figures 3, 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 𝐸_±. See Figures 7 and 8 in Section 4.

(a)

(b)

Figure 3 

(a) Typical points 𝑁1, 𝑁2, 𝑁3, 𝑁4, 𝑁5, and 𝑁6 in the parameter space near the point 𝑅. (b) Phase portraits of system (1.1) for the flow restricted to the center manifold and its continuation. For parameters at 𝑁1 the equilibrium is asymptotically stable; for parameters at 𝑁2 the equilibrium is a weak stable focus (Hopf point with 𝑙1 negative); for parameters at 𝑁3 the equilibrium is unstable and a stable limit cycle appears from a Hopf bifurcation; for parameters at 𝑁4 the equilibrium is an weak unstable focus (Hopf point with 𝑙1 positive) and there is a stable limit cycle; for parameters at 𝑁5 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 𝑁6 the equilibrium is asymptotically stable and the two cycles collide, giving rise to a nondegenerate fold bifurcation of the cycles.

(a)

(b)

Figure 4 

(a) Typical points 𝑀1, 𝑀2, 𝑀3, 𝑀4, 𝑀5, and 𝑀6 in the parameter space near the point 𝑇. (b) Phase portraits of system (1.1) for the flow restricted to the center manifold and its continuation. For parameters at 𝑀1 the equilibrium is unstable; for parameters at 𝑀2 the equilibrium is an weak unstable focus (Hopf point with 𝑙1 positive); for parameters at 𝑀3 the equilibrium is stable and an unstable limit cycle appears from a Hopf bifurcation; for parameters at 𝑀4 the equilibrium is an weak stable focus (Hopf point with 𝑙1 negative) and there is an unstable limit cycle; for parameters at 𝑀5 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 𝑀6 the equilibrium is unstable and the two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles.

We go further and perform a numerical study of the continuation of periodic orbits which arise from the Hopf bifurcations at the points 𝐸_±. 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 Section 2 through a linear analysis of system (1.1) we obtain the Hopf hypersurfaces for 𝐸₀ and 𝐸_±. 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) where

𝐱=(𝑥,𝑦,𝑧)∈ℝ³,𝜻=(𝜇,𝜈,𝛼,𝛽)∈(−∞,∞)×(0,∞)×[1,∞)×(0,∞). (2.2)

The origin 𝐸₀ is an equilibrium point of system (1.1) for all values of the parameters. The characteristic polynomial of the Jacobian matrix of the function 𝑓 given in (2.1) at 𝐸₀ is given by

𝑝(𝜆)=𝜆³+(𝛼−𝜇𝜈)𝜆²+(𝛽−𝜈−𝛼𝜇𝜈)𝜆−𝜇𝛽𝜈. (2.3) If 𝜇>0 then the term −𝜇𝛽𝜈 in the above polynomial is negative and this implies that 𝐸₀ 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 polynomial 𝐿(𝜆)=𝑝₀𝜆³+𝑝₁𝜆²+𝑝₂𝜆+𝑝₃, 𝑝₀>0, with real coefficients, has all roots with negative real parts if and only if the numbers 𝑝₁, 𝑝₂, 𝑝₃ are positive and the inequality 𝑝₁𝑝₂>𝑝₀𝑝₃ is satisfied.

The following proposition is a direct consequence of Lemma 2.1.

Proposition 2.2. Consider 𝜇<0 in system (1.1). If 𝛼𝜇≤−1 then the equilibrium point 𝐸₀ is asymptotically stable for all 𝛽>0 and 𝜈>0. If 𝛼𝜇>−1 and 𝛽>𝛽_1𝑐=𝜈(𝛼−𝜇𝜈)(1+𝛼𝜇)𝛼, (2.4) then system (1.1) has an asymptotically stable equilibrium point at 𝐸₀. If 𝛼𝜇>−1 and 0<𝛽<𝛽_1𝑐 then 𝐸₀ 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 𝐸_±=(±√𝜇,0,±√𝜇) appear. System (1.1) is invariant under the change of coordinates (𝑥,𝑦,𝑧)→(−𝑥,−𝑦,−𝑧). So the stability of the equilibrium 𝐸₋ can be obtained from the stability of 𝐸₊.

The characteristic polynomial of the Jacobian matrix of the function 𝑓 given in (2.1) at 𝐸₊ is given by

𝑝(𝜆)=𝜆³+(𝛼+2𝜇𝜈)𝜆²+(𝛽−𝜈+2𝛼𝜇𝜈)𝜆+2𝜇𝛽𝜈. (2.5) 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 𝐸₊ is asymptotically stable for all 𝛽>0 and 𝜈>0. If 2𝛼𝜇<1 and 𝛽>𝛽_2𝑐=𝜈(𝛼+2𝜇𝜈)(1−2𝛼𝜇)𝛼, (2.6) then system (1.1) has an asymptotically stable equilibrium point at 𝐸₊. If 2𝛼𝜇<1 and 0<𝛽<𝛽_2𝑐 then 𝐸₊ is unstable.

The equations 𝛽=𝛽_1𝑐 and 𝛽=𝛽_2𝑐 in (2.4) and (2.6) give the equations of the Hopf hypersurfaces 𝐻₁ and 𝐻₂ in the parameter space (𝜇,𝜈,𝛼,𝛽). These equations will be used in Section 4 in the study of degenerate Hopf bifurcations which occur at the equilibria 𝐸₀ and 𝐸_± 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 𝑙₁ and 𝑙₂, 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) where 𝐱∈ℝ³ and 𝜻∈ℝ⁴ are, respectively, vectors representing phase variables and control parameters. Assume that 𝑓 is of class 𝐶^∞ in ℝ³×ℝ⁴. Suppose that (3.1) has an equilibrium point 𝐱=𝐱_𝟎 at 𝜻=𝜻_𝟎 and, denoting the variable 𝐱−𝐱_𝟎 also by 𝐱, write

𝐹(𝐱)=𝑓(𝐱,𝜻_𝟎) (3.2) as

1𝐹(𝐱)=𝐴𝐱+21𝐵(𝐱,𝐱)+61𝐶(𝐱,𝐱,𝐱)+124𝐷(𝐱,𝐱,𝐱,𝐱)++1120𝐸(𝐱,𝐱,𝐱,𝐱,𝐱)1720𝐾(𝐱,𝐱,𝐱,𝐱,𝐱,𝐱)+5040𝐿(𝐱,𝐱,𝐱,𝐱,𝐱,𝐱,𝐱)+𝑂(‖𝐱‖⁸), (3.3) where 𝐴=𝑓_𝐱(0,𝜻_𝟎) and, for 𝑖=1,2,3,

𝐵_𝑖(𝐱,𝐲)=³_𝑗,𝑘=1𝜕²𝐹_𝑖(𝜉)𝜕𝜉_𝑗𝜕𝜉_𝑘|_𝜉=0𝑥_𝑗𝑦_𝑘,𝐶_𝑖(𝐱,𝐲,𝐳)=³_{𝑗,𝑘,𝑙=1}𝜕³𝐹_𝑖(𝜉)𝜕𝜉_𝑗𝜕𝜉_𝑘𝜕𝜉_𝑙|_𝜉=0𝑥_𝑗𝑦_𝑘𝑧_𝑙, (3.4) and so on for 𝐷_𝑖, 𝐸_𝑖, 𝐾_𝑖, and 𝐿_𝑖.

Suppose that (𝐱_𝟎,𝜻_𝟎) is an equilibrium point of (3.1) where the Jacobian matrix 𝐴 has a pair of purely imaginary eigenvalues 𝜆_2,3=±𝑖𝜔₀, 𝜔₀>0, and admits no other eigenvalue with zero real part. Let 𝑇^𝑐 be the generalized eigenspace of 𝐴 corresponding to 𝜆_2,3. By this it is meant the largest subspace invariant by 𝐴 on which the eigenvalues are 𝜆_2,3.

Let 𝑝,𝑞∈ℂ³ be vectors such that

𝐴𝑞=𝑖𝜔₀𝑞,𝐴^⊤𝑝=−𝑖𝜔₀𝑝,⟨𝑝,𝑞⟩=³_𝑖=1𝑝_𝑖𝑞_𝑖=1, (3.5) where 𝐴^⊤ is the transpose of the matrix 𝐴. Any vector 𝑦∈𝑇^𝑐 can be represented as 𝑦=𝑤𝑞+𝑤𝑞, where 𝑤⟨𝑝,𝑦⟩∈ℂ. The two-dimensional center manifold associated to the eigenvalues 𝜆_2,3=±𝑖𝜔₀ can be parameterized by the variables 𝑤 and 𝑤 by means of an immersion of the form 𝐱=𝐻(𝑤,𝑤), where 𝐻∶ℂ²→ℝ³ has a Taylor expansion of the form

𝐻(𝑤,𝑤)=𝑤𝑞+𝑤𝑞+_{2≤𝑗+𝑘≤7}1ℎ𝑗!𝑘!^𝑗𝑘𝑤^𝑗𝑤^𝑘+𝑂(|𝑤|⁸), (3.6) with ℎ_𝑗𝑘∈ℂ³ and ℎ_𝑗𝑘=ℎ_𝑘𝑗. Substituting this expression into (3.1) we obtain the following differential equation:

𝐻_𝑤𝑤^′+𝐻_𝑤𝑤^′=𝐹(𝐻(𝑤,𝑤)), (3.7) where 𝐹 is given by (3.2). The complex vectors ℎ_𝑖𝑗 are obtained solving the system of linear equations defined by the coefficients of (3.7), taking into account the coefficients of 𝐹 (see [9, Remark  3.1, page 27]), so that system (3.7), on the chart 𝑤 for a central manifold, writes as follows:

𝑤^′=𝑖𝜔₀1𝑤+2𝐺₂₁𝑤|𝑤|²+1𝐺12³²𝑤|𝑤|⁴+1𝐺144⁴³𝑤|𝑤|⁶+𝑂(|𝑤|⁸), (3.8) with 𝐺_𝑗𝑘∈ℂ.

The first Lyapunov coefficient 𝑙₁ is defined by

𝑙₁=12Re𝐺₂₁, (3.9) where 𝐺₂₁=⟨𝑝,ℋ₂₁⟩, and ℋ₂₁=𝐶(𝑞,𝑞,𝑞)+𝐵(𝑞,ℎ₂₀)+2𝐵(𝑞,ℎ₁₁).

The second Lyapunov coefficient is defined by

𝑙₂=112Re𝐺₃₂, (3.10) where 𝐺₃₂=⟨𝑝,ℋ₃₂⟩. The expression for ℋ₃₂ can be found in [9, equation (36), page 28].

The third Lyapunov coefficient is defined by

𝑙₃=1144Re𝐺₄₃, (3.11) where 𝐺₄₃=⟨𝑝,ℋ₄₃⟩. The expression for ℋ₄₃ can be found in [9, equation (44), page 30].

A Hopf point (𝐱_𝟎,𝜻_𝟎) of system (3.1) is an equilibrium point where the Jacobian matrix 𝐴 has a pair of purely imaginary eigenvalues 𝜆_2,3=±𝑖𝜔₀, 𝜔₀>0, and the other eigenvalue 𝜆₁≠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 called transversal if the parameter dependent complex eigenvalues cross the imaginary axis with nonzero derivative. In a neighborhood of a transversal Hopf point with 𝑙₁≠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 𝑤^′=(𝜂+𝑖𝜔)𝑤+𝑙₁𝑤|𝑤|², where 𝑤∈ℂ, 𝜂, 𝜔, and 𝑙₁ are real functions having derivatives of arbitrary higher order, which are continuations of 0, 𝜔₀, and the first Lyapunov coefficient at the Hopf point. See [8, 9] for details. When 𝑙₁<0(𝑙₁>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.

A Hopf point of codimension 2 is a Hopf point where 𝑙₁ vanishes. It is called transversal if {𝜂=0} and {𝑙₁=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 𝑙₂≠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 𝑤^′=(𝜂+𝑖𝜔₀)𝑤+𝜏𝑤|𝑤|²+𝑙₂𝑤|𝑤|⁴, where 𝜂 and 𝜏 are unfolding parameters. See [8, 9]. The bifurcation diagrams for 𝑙₂≠0 can be found in [8, page 313] and in [16].

A Hopf point of codimension 3 is a Hopf point of codimension 2 where 𝑙₂ vanishes. It is called transversal if {𝜂=0}, {𝑙₁=0}, and {𝑙₂=0} have transversal intersections. In a neighborhood of a transversal Hopf point of codimension 3 with 𝑙₃≠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 𝑤^′=(𝜂+𝑖𝜔₀)𝑤+𝜏𝑤|𝑤|²+𝜈𝑤|𝑤|⁴+𝑙₃𝑤|𝑤|⁶, where 𝜂, 𝜏, and 𝜈 are unfolding parameters. The bifurcation diagram for 𝑙₃≠0 can be found in Takens [16].

4. Hopf Bifurcations in System (1.1)

4.1. Hopf Bifurcations at 𝐸₀

In this subsection we study the stability of 𝐸₀ under the conditions 𝛽=𝛽_1𝑐 and −1<𝛼𝜇<0, that is, on the Hopf hypersurface 𝐻₁ complementary to the range of validity of Proposition 2.2.

The Jacobian matrix of the function 𝑓 given in (2.1) at 𝐸₀ is given by

𝐴=𝐷𝑓(𝐸₀)=𝜇𝜈𝜈01−𝛼−10𝛽_1𝑐0. (4.1) Then, using the notation of the previous section (see expression (3.3)) the multilinear symmetric functions corresponding to 𝑓 can be written as

𝐵(𝐱,𝐲)=𝐷(𝐱,𝐲,𝐳,𝐮)=𝐸(𝐱,𝐲,𝐳,𝐮,𝐯)=(0,0,0),𝐶(𝐱,𝐲,𝐳)=(−6𝜈𝑥₁𝑦₁𝑧₁,0,0). (4.2)

The eigenvalues of 𝐴 are

𝜆₁=−𝛼+𝜇𝜈<0,𝜆_2,3=±𝑖𝜔₀=±𝑖𝜈−𝜇(1+𝛼𝜇)𝛼, (4.3) and from (3.5) one has 𝑞=(−𝑖𝜔₀𝜈𝜇𝜈−𝑖𝜔₀,𝑖𝜔₀,𝛽_1𝑐),1𝑝=𝜌₀(𝑖𝜔₀(𝜇𝜈−𝑖𝜔₀),−𝑖𝜔₀(𝜇𝜈−𝑖𝜔₀)²,(𝜇𝜈−𝑖𝜔₀)²), (4.4) where 𝜌₀=𝛽_1𝑐(𝜇𝜈−𝑖𝜔₀)²+(𝜈+(𝜇𝜈−𝑖𝜔₀)²)𝜔²₀. (4.5)

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 𝐸₀ is given by 𝑙₁(𝜇,𝜈,𝛼,𝛽_1𝑐)=−3𝛼(𝛼𝜇+1)𝜈³(𝛼²−2𝜇𝜈𝛼−𝜈)2(𝛼³−2𝜇𝜈𝛼²−𝜇𝜈²). (4.6) If 𝑠(𝜇,𝜈,𝛼)=−𝛼²+2𝜇𝜈𝛼+𝜈 (4.7) is different from zero then system (1.1) has a transversal Hopf point at 𝐸₀ for 𝛽=𝛽_1𝑐 and −1<𝛼𝜇<0.
More precisely, if (𝜇,𝜈,𝛼,𝛽_1𝑐)∈𝑆∪𝑈 (see Figure 2(a)), then system (1.1) has a Hopf point of codimension 1 at 𝐸₀. If (𝜇,𝜈,𝛼,𝛽_1𝑐)∈𝑆 then 𝐸₀ is asymptotically stable and for each 𝛽<𝛽_1𝑐, but close to 𝛽_1𝑐, there exists a stable periodic orbit near the unstable equilibrium point 𝐸₀. If (𝜇,𝜈,𝛼,𝛽_1𝑐)∈𝑈 then 𝐸₀ is unstable and for each 𝛽>𝛽_1𝑐, but close to 𝛽_1𝑐, there exists an unstable periodic orbit near the asymptotically stable equilibrium point 𝐸₀.

Proof. As the function 𝐵(𝐱,𝐲)=(0,0,0) then the expression 𝐺₂₁ in (3.9) reduces to 𝐺₂₁⟨𝑝,𝐶(𝑞,𝑞,𝑞)⟩. Performing the calculations we get 𝐺₂₁=−6𝜈⁴𝜔⁴₀(𝜇²𝜈²+𝜔²₀)(𝛽_1𝑐(𝜇𝜈−𝑖𝜔₀)²+𝜔²₀(𝜇²𝜈²−2𝑖𝜇𝜔₀𝜈+𝜈−𝜔²₀.)) (4.8) Substituting the expressions of 𝜔₀ and 𝛽_1𝑐 into the above expression of 𝐺₂₁ and taking the real part we have, from (3.9), the expression of the first Lyapunov coefficient 𝑙₁ given in (4.6). Note that the sign of 𝑙₁ is given by the sign of the function 𝑠 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 𝛽=𝛽_1𝑐 verifies 𝜂^′(𝛽_1𝑐)=Re⟨𝑝,𝑑𝐴|𝑑𝛽_{𝛽=𝛽1𝑐}𝛼𝑞⟩=−²2(𝛼³−2𝜇𝜈𝛼²−𝜇𝜈²)<0. (4.9) Therefore, the transversality condition at the Hopf point holds.
The sets 𝑆 and 𝑈 illustrated in Figure 2(a) are defined implicitly as {𝑙₁<0} and {𝑙₁>0}, respectively. The theorem is proved.

Theorem 4.2. Consider the four-parameter family of differential equations (1.1) restricted to 𝛽=𝛽_1𝑐 and −1<𝛼𝜇<0. Then the second Lyapunov coefficient 𝑙₂ associated to the equilibrium 𝐸₀ is negative on the Hopf hypersurface where the first Lyapunov coefficient 𝑙₁ vanishes.

Proof. Following the notation introduced in Section 3, we have 𝐵(𝐱,𝐲)=𝐷(𝐱,𝐲,𝐳,𝐮)=𝐸(𝐱,𝐲,𝐳,𝐮,𝐯)=(0,0,0). Then ℋ₃₂=3𝐶(𝑞,𝑞,ℎ₂₁)+6𝐶(𝑞,𝑞,ℎ₂₁)+3𝐶(𝑞,𝑞,ℎ₃₀). (4.10) The expressions of the complex vectors ℎ₂₁ and ℎ₃₀ 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 coefficient 𝑙₂ associated to 𝐸₀ can be written as

𝑙₂(𝜇,𝜈,𝛼,𝛽_1𝑐)=−3𝛼²(𝛼𝜇+1)²𝜈⁵𝑁(𝜇,𝜈,𝛼,𝛽_1𝑐)2𝜇(−𝛼³+2𝜇𝜈𝛼²+𝜇𝜈²)³(−𝛼³+2𝜇𝜈𝛼²+8𝜇²𝜈²𝛼+9𝜇𝜈²), (4.11) where 𝑁(𝜇,𝜈,𝛼,𝛽_1𝑐) has the form

𝛼¹⁰−8𝜇𝜈𝛼⁹+𝜈(16𝜇²𝜈−1)𝛼⁸+2𝜇𝜈²(8𝜇²𝜈−7)𝛼⁷+2𝜇²𝜈³(61−40𝜇²𝜈)𝛼⁶+2𝜇𝜈³(32𝜈²𝜇⁴−112𝜈𝜇²+13)𝛼⁵−𝜇²𝜈⁴(64𝜈𝜇²+19)𝛼⁴+40𝜇³𝜈⁵(10𝜇²𝜈−9)𝛼³+𝜇²𝜈⁵(−192𝜈²𝜇⁴+724𝜈𝜇²−121)𝛼²−8𝜇³𝜈⁶(43𝜇²𝜈−34)𝛼−142𝜇⁴𝜈⁷. (4.12) Along the surface {𝑙₁=0} we have, from (4.7) of Theorem 4.1, 𝑠(𝜇,𝜈,𝛼)=−𝛼²+2𝜇𝜈𝛼+𝜈=0, or equivalently, 𝜈𝛼²/(1+2𝛼𝜇). Substituting this expression into (4.11) one has

𝑙₂|_{𝑙1=0}=−9𝛼¹¹(𝛼𝜇+1)(4𝛼𝜇−1)(2𝛼𝜇+1)⁴(8𝛼𝜇−1)<0. (4.13)

Here the transversality condition is equivalent to the transversality of the hypersurfaces 𝛽=𝛽_1𝑐 and 𝑙₁(𝜇,𝜈,𝛼,𝛽)=0. The bifurcation diagram for a typical point 𝑅 on the surface {𝑙₁=0} is depicted in Figure 3.

4.2. Hopf Bifurcations at 𝐸_±

In this subsection we study the stability of 𝐸_± under the conditions 𝛽=𝛽_2𝑐 and 0<2𝛼𝜇<1, that is, on the Hopf hypersurface 𝐻₂ complementary to the range of validity of Proposition 2.3.

The Jacobian matrix of the function 𝑓 given in (2.1) calculated at 𝐸₊ is given by

𝐴=𝐷𝑓(𝐸₊)=−2𝜇𝜈𝜈01−𝛼−10𝛽_2𝑐0. (4.14) Then, using the notation of Section 3 (see expression (3.3)) the multilinear symmetric functions corresponding to 𝑓 can be written as

𝐵(𝐱,𝐲)=(−6𝜈√𝜇𝑥₁𝑦₁,0,0),𝐶(𝐱,𝐲,𝐳)=(−6𝜈𝑥₁𝑦₁𝑧₁,0,0),𝐷=𝐸=𝐾=𝐿=(0,0,0). (4.15)

The eigenvalues of 𝐴 are

𝜆₁=−𝛼−2𝜇𝜈<0,𝜆_2,3=±𝑖𝜔₀=±𝑖𝜈2𝜇(1−2𝛼𝜇)𝛼, (4.16) and from (3.5) one has

𝑞=(𝑖𝜔₀𝜈2𝜇𝜈+𝑖𝜔₀,𝑖𝜔₀,𝛽_2𝑐),1𝑝=𝜌₀(−𝑖𝜔₀(2𝜇𝜈+𝑖𝜔₀),−𝑖𝜔₀(2𝜇𝜈+𝑖𝜔₀)²,(2𝜇𝜈+𝑖𝜔₀)²), (4.17)(4.18) where

𝜌₀=𝛽_2𝑐(2𝜇𝜈+𝑖𝜔₀)²+(𝜈+(2𝜇𝜈+𝑖𝜔₀)²)𝜔²₀. (4.19)

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 𝐸₊ is given by 𝑙₁(𝜇,𝜈,𝛼,𝛽_2𝑐)=−3𝛼(2𝛼𝜇−1)𝜈³𝑔(𝜇,𝜈,𝛼)(𝛼³+4𝜇𝜈𝛼²+2𝜇𝜈²)(𝛼³+4𝜇𝜈𝛼²−12𝜇²𝜈²𝛼+8𝜇𝜈²), (4.20) where 𝑔(𝜇,𝜈,𝛼)=𝛼⁵+8𝜇𝜈𝛼⁴+𝜈(4𝜇²𝜈−1)𝛼³+2𝜇𝜈²(5−24𝜇²𝜈)𝛼²+80𝜇²𝜈³𝛼−20𝜇𝜈³. (4.21) If 𝑔(𝜇,𝜈,𝛼) is different from zero then system (1.1) has a transversal Hopf point at 𝐸₊ for 𝛽=𝛽_2𝑐 and 0<2𝛼𝜇<1.
More precisely, if (𝜇,𝜈,𝛼,𝛽_2𝑐)∈𝑆∪𝑈 (see Figure 2(b)) then system (1.1) has a Hopf point of codimension 1 at 𝐸₊. That is, if (𝜇,𝜈,𝛼,𝛽_2𝑐)∈𝑆 then 𝐸₊ is asymptotically stable and for each 𝛽<𝛽_2𝑐, but close to 𝛽_2𝑐, there exists a stable periodic orbit near the unstable equilibrium point 𝐸₊; if (𝜇,𝜈,𝛼,𝛽_2𝑐)∈𝑈 then 𝐸₊ is unstable and for each 𝛽>𝛽_2𝑐, but close to 𝛽_2𝑐, there exists an unstable periodic orbit near the asymptotically stable equilibrium point 𝐸₊.

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 𝑞 in (4.17), the vector 𝑝 in (4.18) and the multilinear functions 𝐵 and 𝐶. 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 −𝑔 defined in (4.21) since the denominator of 𝑙₁ and the expression 3(2𝛼𝜇−1)𝜈³ are negative. In Figure 2(b) the surface {𝑙₁=0} and the regions 𝑆 and 𝑈 are illustrated. These regions are defined implicitly as {𝑙₁<0} and {𝑙₁>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 𝛽=𝛽_2𝑐 verifies 𝜂^′(𝛽_2𝑐)=Re⟨𝑝,𝑑𝐴|𝑑𝛽_{𝛽=𝛽2𝑐}𝛼𝑞⟩=−²2(𝛼³+4𝜇𝜈𝛼²+2𝜇𝜈²)<0. (4.22) 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 𝛽=𝛽_2𝑐 and 0<2𝛼𝜇<1. Then the second Lyapunov coefficient 𝑙₂ associated to the equilibrium 𝐸₊ is given by 𝑙₂(𝜇,𝜈,𝛼,𝛽_2𝑐𝛼)=²(2𝛼𝜇−1)𝜈⁵ℎ(𝜇,𝜈,𝛼)4𝜇(𝛼³+4𝜇𝜈𝛼²−32𝜇²𝜈²𝛼+18𝜇𝜈²,)𝑑(𝜇,𝜈,𝛼) (4.23) where the function ℎ(𝜇,𝜈,𝛼) is given by 21𝜇𝛼²⁰+(588𝜇²𝜈−4)𝛼¹⁹+3𝜇𝜈(1876𝜇²𝜈−53)𝛼¹⁸+4𝜈(3192𝜈²𝜇⁴−294𝜈𝜇²+1)𝛼¹⁷+𝜇𝜈²(−121296𝜈²𝜇⁴+14156𝜈𝜇²−65)𝛼¹⁶−12𝜇²𝜈³(59248𝜈²𝜇⁴−18704𝜈𝜇²+473)𝛼¹⁵+𝜇𝜈³(114240𝜈³𝜇⁶+716976𝜈²𝜇⁴−62932𝜈𝜇²+199)𝛼¹⁴+16𝜇²𝜈⁴(477120𝜈³𝜇⁶−165272𝜈²𝜇⁴+769𝜈𝜇²+119)𝛼¹³+4𝜇³𝜈⁵(1634304𝜈³𝜇⁶−4121136𝜈²𝜇⁴+751264𝜈𝜇²−16945)𝛼¹²−16𝜇²𝜈⁵(2225664𝜈⁴𝜇⁸−327552𝜈³𝜇⁶−585624𝜈²𝜇⁴+53587𝜈𝜇²−262)𝛼¹¹−4𝜇³𝜈⁶(6967296𝜈⁴𝜇⁸−29565184𝜈³𝜇⁶+7648944𝜈²𝜇⁴+81592𝜈𝜇²−11129)𝛼¹⁰+16𝜇⁴𝜈⁷(4644864𝜈⁴𝜇⁸+773376𝜈³𝜇⁶−8268704𝜈²𝜇⁴+1627004𝜈𝜇²−41405)𝛼⁹+4𝜇³𝜈⁷(−80474112𝜈⁴𝜇⁸[phone]𝜈³𝜇⁶[phone]𝜈²𝜇⁴−1903712𝜈𝜇²+17531)𝛼⁸+16𝜇⁴𝜈⁸(5308416𝜈⁴𝜇⁸[phone]𝜈³𝜇⁶−13584144𝜈²𝜇⁴−51152𝜈𝜇²+39713)𝛼⁷+32𝜇⁵𝜈⁹(1327104𝜈⁴𝜇⁸−10505088𝜈³𝜇⁶−13899000𝜈²𝜇⁴+4382138𝜈𝜇²−130059)𝛼⁶−16𝜇⁴𝜈⁹(11446272𝜈⁴𝜇⁸−35076672𝜈³𝜇⁶−10290192𝜈²𝜇⁴+2504248𝜈𝜇²−37769)𝛼⁵+32𝜇⁵𝜈¹⁰(10175616𝜈³𝜇⁶−15588272𝜈²𝜇⁴−317442𝜈𝜇²+129593)𝛼⁴−64𝜇⁶𝜈¹¹(4720704𝜈²𝜇⁴−3846404𝜈𝜇²+142349)𝛼³+128𝜇⁵𝜈¹¹(1199432𝜈²𝜇⁴−495365𝜈𝜇²+13321)𝛼²−256𝜇⁶𝜈¹²(157944𝜇²𝜈−25841)𝛼+4308992𝜇⁷𝜈¹³, (4.24) and the function 𝑑(𝜇,𝜈,𝛼) in the denominator is given by (𝛼⁶+8𝜇𝜈𝛼⁵+4𝜇²𝜈²𝛼⁴+2𝜇𝜈²(5−24𝜇²𝜈)𝛼³+40𝜇²𝜈³𝛼²−24𝜇³𝜈⁴𝛼+16𝜇²𝜈⁴)³. (4.25) If (𝜇,𝜈,𝛼,𝛽_2𝑐)∈𝐶₁∪𝐶₂ (see Figure 2(b)), then system (1.1) has a transversal Hopf point of codimension 2 at 𝐸₊. More specifically, if (𝜇,𝜈,𝛼,𝛽_2𝑐)∈𝐶₁ then the Hopf point at 𝐸₊ is asymptotically stable and the bifurcation diagram is illustrated in Figure 3 for a typical point 𝑅. If (𝜇,𝜈,𝛼,𝛽_2𝑐)∈𝐶₂ then the Hopf point at 𝐸₊ is unstable and the bifurcation diagram is drawn in Figure 4 for a typical point 𝑇.

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 𝑞 in (4.17), the vector 𝑝 in (4.18) and the multilinear functions 𝐵, 𝐶, 𝐷, and 𝐸. The calculations in this proof, corroborated by Computer Algebra, have been posted in [17].
In Figure 2 we present a geometric synthesis interpreting the long calculations involved in this proof. The sign of −ℎ(𝜇,𝜈,𝛼) gives the sign of the second Lyapunov coefficient since the function 𝑑(𝛼,𝜇,𝜈) in the denominator of 𝑙₂ is positive. The graph of ℎ(𝜇,𝜈,𝛼)=0 on the surface {𝑙₁=0}, which determines the curve 𝐶, and the signs of the first and second Lyapunov coefficients are also illustrated. As follows, 𝑙₂<0 on the open region on the surface {𝑙₁=0}, denoted by 𝐶₁. On this region a typical reference point 𝑅 is depicted. Also 𝑙₂>0 on the open region on the surface {𝑙₁=0}, denoted by 𝐶₂. This region contains the typical reference point, denoted by 𝑇.
The bifurcation diagrams of system (1.1) at the points 𝑅 and 𝑇 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 curve 𝐶 given by the intersection of the surfaces {𝑙₁=0} and {𝑙₂=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 𝐶. So the bifurcation diagrams depicted in Figures 7 and 8 are essentially the same for all values of 𝛼≥1.

Figure 5 

Graph of 𝑙3 with respect to 𝛼.

All the calculations presented here are illustrated in Figure 6 for the cases 𝛼=2 (a) and 𝛼=3 (b). In this figure are depicted the surfaces 𝛽=𝛽_2𝑐, the solid curves where 𝑙₁=0, the regions 𝑆 and 𝑈 where 𝑙₁<0 and 𝑙₁>0, respectively, the dashed curves where 𝑙₂=0 and the arcs 𝐶₁ and 𝐶₂ where 𝑙₂<0 and 𝑙₂>0, respectively. The bifurcation diagrams for typical points 𝑅₁ and 𝑅₂ are depicted in Figure 3 while the bifurcation diagrams for typical points 𝑇₁ and 𝑇₂ are shown in Figure 4. In the next theorem we analyze the sign of the third Lyapunov coefficient at the point 𝑄₁, that is, for the case 𝛼=2 or equivalently for 𝑅=𝑅_𝑝. For the case 𝛼=3 the calculations are very similar.

(a)

(b)

Figure 6 

(a) Hopf surface 𝛽=𝛽𝑐 and curves {𝑙1=0} (solid) and {𝑙2=0} (dashed) for 𝛼=2. (b) Hopf surface 𝛽=𝛽𝑐 and curves {𝑙1=0} (solid) and {𝑙2=0} (dashed) for 𝛼=3.

Figure 7 

Bifurcation diagram of system (1.1) for 𝛼=2 at the point 𝑄1.

(a)

(b)

Figure 8 

(a) Typical points 𝑉1, 𝑉2, 𝑉3, 𝑉4, 𝑉5, 𝑉6, 𝑉7, 𝑉8, and 𝑉9 in the parameter space near the point 𝑅1 (see Figure 7). (b) Phase portraits of system (1.1) for the flow restricted to the center manifold and its continuation at the points 𝑉𝑖, 𝑖=1,…,9 of Figure 8(a).

Theorem 4.6. Consider 𝛼=2 in system (1.1). Then there is a unique point 𝑄₁=𝑄_𝛼=2=(𝜇,𝜈,𝛽_2𝑐) where the surfaces {𝑙₁=0} and {𝑙₂=0} intersect transversally on the Hopf hypersurface 𝐻₂. For the parameter values at the point 𝑄₁, system (1.1) has a transversal Hopf point of codimension 3 at 𝐸₊ which is unstable since 𝑙₃(𝑄₁)>0. The bifurcation diagram of system (1.1) at the point 𝑄₁ is illustrated in Figure 7. In Figure 8 the bifurcation diagram of system (1.1) at a typical point 𝑅₁ of Figure 7 is drawn.

Computer Assisted Proof
The point 𝑄₁ is the intersection of the surfaces {𝑙₁=0} and {𝑙₂=0} on the Hopf hypersurface 𝛽=𝛽_2𝑐 taking into account 𝛼=2. It is defined and obtained as the solution of the equations 𝑔(𝜇,𝜈,𝛼)=0 (see Theorem 4.3) and ℎ(𝜇,𝜈,𝛼)=0 (see Theorem 4.4). The point 𝑄₁ has coordinates 𝛽𝜇=0.01528952519719747820…,𝜈=4.03634524181150646828…,_2𝑐=4.02335387131139691747…. (4.26) The existence and uniqueness of 𝑄₁ 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 point 𝑄₁ take five decimal round-off coordinates 𝜇=0.01528, 𝜈=4.03634, and 𝛽_2𝑐=4.02335 for simplicity. For these values of the parameters one has ℎ𝑝=(−0.10542+0.46290𝑖,0.21084+0.10811𝑖,0.22356−0.43599𝑖),𝑞=(3.78949+0.96718𝑖,0.48359𝑖,4.02335),₁₁ℎ=(−371.10203,0,−371.10203),₂₀ℎ=(82.25866+60.23037𝑖,−1.95710+26.99100𝑖,112.27829+8.14125𝑖),₃₀=,𝐺1209.86465+3423.09545𝑖−366.52849+1469.19120𝑖4074.40090+1016.46677𝑖₂₁ℎ=−2878.05753𝑖,₂₁=,ℎ13925.76364+2288.71956𝑖−623.13603−1297.11421𝑖13152.94476+5184.29327𝑖₄₀=,ℎ−83816.60065+166465.53895𝑖−47712.54683+79760.69093𝑖165895.88420+99238.29710𝑖₃₁=2.26342⋅10⁶+1.12088⋅10⁶𝑖−87119.52533+376583.80975𝑖2.56884⋅10⁶,ℎ+435080.93606𝑖₂₂=(−8.42248⋅10⁶,0,−8.42248⋅10⁶𝐺),₃₂=−9.92235⋅10⁷ℎ𝑖,₃₂=5.84180⋅10⁸+1.09060⋅10⁸𝑖−1.70175⋅10⁷−6.11162⋅10⁷𝑖5.51875⋅10⁸+2.34141⋅10⁸𝑖,ℎ₄₁=1.25475⋅10⁸+1.68063⋅10⁸𝑖−1.64672⋅10⁷+6.39292⋅10⁷𝑖2.25787⋅10⁸+5.77662⋅10⁷𝑖,ℎ₄₂=1.93748⋅10¹¹+9.38389⋅10¹⁰𝑖−7.94067⋅10⁹+2.60031⋅10¹⁰𝑖2.15396⋅10¹¹+4.67300⋅10¹⁰𝑖,ℎ₃₃=(−7.28682⋅10¹¹,0,−7.28682⋅10¹¹𝐺),₄₃=1.17553⋅10¹⁰−1.14289⋅10¹³𝑖. (4.27) From (3.9), (3.10), (3.11), and (4.27) one has 𝑙₁(𝑄₁)=0,𝑙₂(𝑄₁)=0,𝑙₃(𝑄₁1)=144Re𝐺₄₃=1144(1.17553⋅10¹⁰). (4.28) 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 functions 𝑙₁, given in (4.20), and 𝑙₂, given in (4.23), at the point 𝑄₁ are, respectively, (181.95647,−14.19675),(−5.26376⋅10⁶,4.12653⋅10⁶). (4.29) The transversality condition at 𝑄₁ 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⋅10⁸>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.
Define 𝑙₀=𝛽−𝛽_2𝑐. The open region in the parameter space where system (1.1) has three small periodic orbits bifurcating from the equilibria 𝐸_± can be described by 𝑙₂<0, 𝑙₁>0, and 𝑙₀>0 with |𝑙₂|≫𝑙₁≫𝑙₀>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 𝑉₁ the equilibrium is unstable; for parameters at 𝑉₂ the equilibrium is an weak unstable focus (Hopf point with positive 𝑙₁); for parameters at 𝑉₃ the equilibrium is stable and an unstable limit cycle appears from the Hopf bifurcation; for parameters at 𝑉₄ the equilibrium is an weak stable focus (Hopf point with negative 𝑙₁) and there is an unstable limit cycle; for parameters at 𝑉₅ 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 𝑉₆ the equilibrium is unstable and two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles; for parameters at 𝑉₇ the equilibrium is stable and two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles; for parameters at 𝑉₈ the equilibrium is stable and two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles; for parameters at 𝑉₉ 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 𝐸₀ and 𝐸_±, 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 𝐸_±. 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 at 𝐸₀

For 𝜇<0, system (1.1) has the origin 𝐸₀=(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 𝛽_1𝑐=𝜈(𝛼−𝜇𝜈)(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 𝛽_1𝑐=2.0625. One can observe that when the parameter 𝛽 tends to the critical value 𝛽=𝛽_1𝑐 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 𝛽=𝛽_1𝑐 is very small (observe the range of the state variables (𝑥,𝑦,𝑧) in Figure 10). The numerical analysis performed suggests that such small limit cycle which birth in the Hopf bifurcation for 𝛽=𝛽_1𝑐 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 𝛽_1𝑐 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].

(a)

(b)

Figure 9 

Phase portrait of system (1.1) near the equilibrium point 𝐸0 for the parameter values 𝜈=3, 𝜇=−0.25, 𝛼=2, and (a) 𝛽=𝛽1𝑐+0.5=2.5625; (b) 𝛽=𝛽1𝑐+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 𝐸0 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 𝛽1𝑐=2.0625. Eigenvalues at 𝐸0 are (a) 𝜆1=−2.624, 𝜆2,3=−0.062±0.853𝑖; (b) 𝜆1=−2.737, 𝜆2,3=−0.006±0.760𝑖.

(a)

(b)

Figure 10 

Limit cycle created in the Hopf bifurcation for system (1.1) at the point 𝐸0. The parameter values are 𝜈=3, 𝜇=−0.25, 𝛼=2, and 𝛽=𝛽1𝑐−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 𝐸0 is an unstable focus and a limit cycle arises around it. Eigenvalues at 𝐸0 are 𝜆1=−2.774,𝜆2,3=0.012±0.728𝑖.

(a)

(b)

Figure 11 

Limit cycle created in the Hopf bifurcation for system (1.1) at the point 𝐸0. The parameter values are 𝜈=3, 𝜇=−0.25, 𝛼=2, and 𝛽=𝛽1𝑐−1=1.0625 in (a); 𝛽=𝛽1𝑐−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 𝛽1𝑐 the solutions tend faster to the limit cycle. Eigenvalues at 𝐸0 are (a) 𝜆1=−2.985, 𝜆2,3=0.117±0.502𝑖; (b) 𝜆1=−3.203, 𝜆2=0.418, 𝜆3=0.034.

5.2. Bifurcations at 𝐸_±

For 𝜇>0, system (1.1) has the origin 𝐸₀ and also the symmetric equilibria 𝐸_±=(±√𝜇,0,±√𝜇). According to Theorem 4.3, the system undergoes a Hopf bifurcation when the parameter 𝛽 crosses the critical value 𝛽=𝛽_2𝑐=𝜈(𝛼+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 𝛽_2𝑐=12. In the figures are shown the 1-dimensional unstable manifolds of 𝐸₀, which spiral to the equilibria 𝐸_± for 𝛽>𝛽_2𝑐. Observe that when the parameter 𝛽 tends to the critical value 𝛽=𝛽_2𝑐 the spiralling behavior of the solutions becomes slower (Figure 12(b)); for 𝛽<𝛽_2𝑐 the unstable manifolds spiral to the symmetric limit cycles arise in the Hopf bifurcation for 𝛽=𝛽_2𝑐 (Figure 13).

(a)

(b)

Figure 12 

Phase portraits of system (1.1) before the Hopf bifurcation at the points 𝐸±. The parameter values are 𝜈=10, 𝜇=0.1, 𝛼=2, and 𝛽=𝛽2𝑐+0.8=12.8 in (a); 𝛽=𝛽2𝑐+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 𝐸0 are (a) 𝜆1=1.951, 𝜆2,3=−1.475±2.093𝑖; (b) 𝜆1=1.987, 𝜆2,3=−1.493±1.976𝑖. Eigenvalues at 𝐸± are (a) 𝜆1=−3.928, 𝜆2,3=−0.035±2.552𝑖; (b) 𝜆1=−3.981, 𝜆2,3=−0.009±2.475𝑖.

(a)

(b)

Figure 13 

Phase portraits of system (1.1) after the Hopf bifurcation at the points 𝐸±. The parameter values are 𝜈=10, 𝜇=0.1, 𝛼=2, and 𝛽=𝛽2𝑐−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 𝐸0 are 𝜆1=2.012, 𝜆2,3=−1.506±1.895𝑖. Eigenvalues at 𝐸± are 𝜆1=−4.018, 𝜆2,3=0.009±2.423𝑖.

We observe that when the parameter 𝛽 moves away from the critical value 𝛽_2𝑐, 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 𝐸₀ 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 𝐸_± for the critical value of parameter 𝛽=𝛽_2𝑐. For the case 𝛼=1 another mechanism giving rise to chaotic attractor is described in [5].

(a)

(b)

Figure 14 

(a) Period doubling bifurcation for the limit cycle created in the Hopf bifurcation at the point 𝐸±. The parameter values are 𝜈=10, 𝜇=0.1, 𝛼=2, and 𝛽=𝛽2𝑐−1.13=10.87. (b) Double symmetric homoclinic loop to the origin. Parameter values 𝜈=10, 𝜇=0.1, 𝛼=2, and 𝛽=𝛽2𝑐−1.214322=10.785678. Initial conditions are (0.0001,0,0.0001) and (−0.0001,0,−0.0001). Eigenvalues at 𝐸0 are (a) 𝜆1=2.073, 𝜆2,3=−1.536±1.697𝑖; (b) 𝜆1=2.079, 𝜆2,3=−1.539±1.678𝑖. Eigenvalues at 𝐸± are (a) 𝜆1=−4.104, 𝜆2,3=0.052±2.300𝑖; (b) 𝜆1=−4.111, 𝜆2,3=0.055±2.289𝑖.

Figure 15 

Strange attractor created after the existence of the double homoclinic loop shown in Figure 14. Parameter values are 𝜈=10, 𝜇=0.1, 𝛼=2, and 𝛽=𝛽2𝑐−1.743=10.257. Initial conditions are (0.0001,0,0.0001) and (−0.0001,0,−0.0001). Eigenvalues at 𝐸0 are 𝜆1=2.115, 𝜆2,3=−1.557±1.556𝑖; eigenvalues at 𝐸± are 𝜆1=−4.161, 𝜆2,3=0.080±2.218𝑖.

5.3. Bifurcation at the Critical Point 𝑄₁ in the Parameter Space

As stated in Theorem 4.6, the most complex Hopf bifurcation at the points 𝐸_± occurs at the point 𝑄₁=(𝜇,𝜈,𝛽_2𝑐) 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 𝑄₁=(𝜇_𝑄1,𝜈_𝑄1,𝛽_𝑄1)=(0.015289525,4.036345241,4.023353871). In Figure 16 we show an approximation of the unstable manifolds of the equilibrium 𝐸₀. These manifolds spiral to the equilibria 𝐸_± if 𝛽>𝛽_𝑄1=4.023353871 and this spiralling behavior become slower as 𝛽 tends to 𝛽_𝑄1=𝛽_2𝑐 from above (for 𝜇=𝜇_𝑄1 and 𝜈=𝜈_𝑄1 fixed).

(a)

(b)

Figure 16 

Phase portraits of system (1.1) before the Hopf bifurcation at the points 𝐸±. The parameter values are near the point 𝑄1 of Theorem 4.6 in the parameter space: 𝜇=𝜇𝑄1=0.015289525, 𝜈=𝜈𝑄1=4.036345241, and 𝛽=𝛽𝑄+0.5=4.523353871 in (a); 𝛽=𝛽𝑄+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 𝐸0 are (a) 𝜆1=0.280, 𝜆2=−1.602, 𝜆3=−0.617; (b) 𝜆1=0.357, 𝜆2=−1.937, 𝜆3=−0.358. Eigenvalues at 𝐸± are (a) 𝜆1=−1.891, 𝜆2,3=−0.115±0.530𝑖; (b) 𝜆1=−2.119, 𝜆2,3=−0.002±0.484𝑖.

For (𝜇,𝛼,𝛽) at the point 𝑄₁ (see Theorem 4.6) the system presents a highly degenerate behavior. In fact, due to the vanishing of the first and second Lyapunov coefficients, 𝑙₁ and 𝑙₂, the solutions near the equilibria 𝐸_± 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 𝐸_± in Theorem 4.6, since 𝑙₃(𝑄₁)>0, the solutions move away from these points very slowly, with repeller rate as 10⁻⁹. Then, it is not detected numerically and the solutions near the equilibria 𝐸_± seem to be periodic as shown in Figure 17. In addition there is also a “big" periodic orbit encircling the equilibria 𝐸_±. It would be very interesting to observe how the real electronic circuit behave for the parameter values at the bifurcation point 𝑄₁.

Figure 17 

Phase portrait of system (1.1) for the parameters at the point 𝑄1=(𝜇𝑄1,𝜈𝑄1,𝛽𝑄1) of Theorem 4.6 in the parameter space: (𝜇,𝜈.𝛽). The solutions near the equilibria 𝐸± 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 𝐸±. Eigenvalues at 𝐸0 are 𝜆1=0.359,𝜆2=−1.942,𝜆3=−0.355. Eigenvalues at 𝐸± are 𝜆1=−2.123, 𝜆2,3=−0.312×10−10±0.483𝑖≈±0.483𝑖.

For 𝛽<𝛽_𝑄1, 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.

(a)

(b)

Figure 18 

Phase portrait of system (1.1) for 𝜇=𝜇𝑄1, 𝜈=𝜈𝑄1, and 𝛽=3.973353871<𝛽𝑄1. The three equilibria are unstable and there exists a “big” stable periodic orbit encircling them. Eigenvalues at 𝐸0 are 𝜆1=0.368, 𝜆2=−1.969, 𝜆3=−0.337. Eigenvalues at 𝐸± are 𝜆1=−2.144, 𝜆2,3=0.010±0.478𝑖.

6. Concluding Remarks

In this paper we analyze the Lyapunov stability of the equilibria 𝐸₀ and 𝐸_± 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 𝐸₀ (resp., 𝐸_±) 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. Figures 3, 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 𝐸_±.

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 Section 6. 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.

Acknowledgments

The first author is supported by CAPES. The second author developed this work under the Project CNPq 473747/2006-5. The third author is partially supported by CNPq under the Project 478544/2007-3. The authors also would like to thank the financial support from PROINTER/PROPe/UNESP and FUNDUNESP.