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

1 Instituto de Sistemas Elétricos e Energia, Universidade Federal de Itajubá, Avenida BPS 1303, Pinheirinho, 37500-903 Itajubá, MG, Brazil 2 Instituto de Ciências Exatas, Universidade Federal de Itajubá, Avenida BPS 1303, Pinheirinho, 37500-903 Itajubá, MG, Brazil 3 Departamento de Matemática, Estatı́stica e Computação, Faculdade de Ciências e Tecnologia, UNESP—Univ Estadual Paulista, Cx. Postal 266, 19060-900 Presidente Prudente, SP, Brazil


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: where x, y, z ∈ R 3 are the state variables and μ ∈ R, α ≥ 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 by the following changes in variables, parameters, and rescaling in time:

1.3
In 1.2 v 1 and v 2 are the voltages across the capacitors C 1 and C 2 , respectively, i N is the current through the nonlinear diode N, i L is the current through the inductor, c 1 and c 2 are the capacitances, L is the inductance, R is the linear resistor, R p is the parallel resistor, α 1 R/R p , and the prime denotes derivatives with respect to the time t.It is assumed that the current through the nonlinear diode is given by i N i N v 1 av 1 bv 3  1 , with a < 0 and b > 0. See Figure 1 and 1 for more details.
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 E 0 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 E 0 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 E 0 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 equilibrium E 0 0, 0, 0 the Hopf hypersurface H 1 is obtained in the space of parameters μ, ν, α, β and the first Lyapunov coefficient l 1 is calculated.It is shown that this coefficient vanishes on a 2-dimensional surface contained in H 1 , giving rise to codimension two Hopf bifurcations see Figure 2 a .Then the second Lyapunov coefficient l 2 is calculated and it is established that this coefficient is always negative on the surface {l 1 0}.
b For the symmetric equilibria E ± ± √ μ, 0, ± √ μ the Hopf hypersurface H 2 is obtained in the space of parameters μ, ν, α, β , the first Lyapunov coefficient l 1 is calculated, and it is shown that this coefficient vanishes on a 2-dimensional surface contained in H 2 see Figure 2 b , giving rise to codimension two bifurcations.The second Lyapunov coefficient l 2 is calculated and it is established that this coefficient also vanishes on a 2-dimensional surface contained in H 2 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 {l 1 0} and {l 2 0} are obtained and found to be positive.
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 E ± .See Figures 7 and 8 in Section 4.
We go further and perform a numerical study of the continuation of periodic orbits which arise from the Hopf bifurcations at the points E ± .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 E 0 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.

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 where The origin E 0 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 E 0 is given by If μ > 0 then the term −μβν in the above polynomial is negative and this implies that E 0 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 L λ p 0 λ 3 p 1 λ 2 p 2 λ p 3 , p 0 > 0, with real coefficients, has all roots with negative real parts if and only if the numbers p 1 , p 2 , p 3 are positive and the inequality p 1 p 2 > p 0 p 3 is satisfied.
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 function f given in 2.1 at E is given by The following proposition is also a direct consequence of Lemma 2.1.
The equations β β 1c and β β 2c in 2.4 and 2.6 give the equations of the Hopf hypersurfaces H 1 and H 2 in the parameter space μ, ν, α, β .These equations will be used in Section 4 in the study of degenerate Hopf bifurcations which occur at the equilibria E 0 and E ± of system 1.1 .

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 l 1 and l 2 , 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 and so on for D i , E i , K i , and L i .Suppose that x 0 , ζ 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 T c 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 .
Let p, q ∈ C 3 be vectors such that Aq iω 0 q, A p −iω 0 p, p, q where A is the transpose of the matrix A. Any vector y ∈ T c can be represented as y wq w q, where w p, y ∈ C. 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 : p, H 43 .The expression for H 43 can be found in 9, equation 44 , page 30 .A Hopf point x 0 , ζ 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 called transversal if the parameter dependent complex eigenvalues cross the imaginary axis with nonzero derivative.In a neighborhood of a transversal Hopf point with l 1 / 0 the dynamic behavior of the system 3.1 , reduced to the family of parameterdependent continuations of the center manifold, is orbitally topologically equivalent to the following complex normal form w η iω w l 1 w|w| 2 , where w ∈ C, η, ω, and l 1 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 l 1 < 0 l 1 > 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 l 1 vanishes.It is called transversal if {η 0} and {l 1 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 l 2 / 0 the dynamic behavior of the system 3.1 reduced to the family of parameterdependent continuations of the center manifold is orbitally topologically equivalent to w η iω 0 w τw|w| 2 l 2 w|w| 4 , where η and τ are unfolding parameters.See 8, 9 .The bifurcation diagrams for l 2 / 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 l 2 vanishes.It is called transversal if {η 0}, {l 1 0}, and {l 2 0} have transversal intersections.In a Mathematical Problems in Engineering neighborhood of a transversal Hopf point of codimension 3 with l 3 / 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 l 3 w|w| 6 , where η, τ, and ν are unfolding parameters.The bifurcation diagram for l 3 / 0 can be found in Takens 16 .

Hopf Bifurcations at E 0
In this subsection we study the stability of E 0 under the conditions β β 1c and −1 < αμ < 0, that is, on the Hopf hypersurface H 1 complementary to the range of validity of Proposition 2.2.
The Jacobian matrix of the function f given in 2.1 at E 0 is given by 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
The eigenvalues of A are and from 3.5 one has where Using these calculations we prove the next two theorems, from which statement a of the introduction follows.
Proof.As the function B x, y 0, 0, 0 then the expression G 21 in 3.9 reduces to G 21 p, C q, q, q .Performing the calculations we get Substituting the expressions of ω 0 and β 1c into the above expression of G 21 and taking the real part we have, from 3.9 , the expression of the first Lyapunov coefficient l 1 given in 4.6 .Note that the sign of l 1 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 Therefore, the transversality condition at the Hopf point holds.
The sets S and U illustrated in Figure 2 a are defined implicitly as {l 1 < 0} and {l 1 > 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 l 2 associated to the equilibrium E 0 is negative on the Hopf hypersurface where the first Lyapunov coefficient l 1 vanishes.

4.10
The expressions of the complex vectors h 21 and h 30 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 l 2 associated to E 0 can be written as where N μ, ν, α, β 1c has the form

Hopf Bifurcations at E ±
In this subsection we study the stability of E ± under the conditions β β 2c and 0 < 2αμ < 1, that is, on the Hopf hypersurface H 2 complementary to the range of validity of Proposition 2.3.
The Jacobian matrix of the function f given in 2.1 calculated at E is given by , and N 6 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 N 1 the equilibrium is asymptotically stable; for parameters at N 2 the equilibrium is a weak stable focus Hopf point with l 1 negative ; for parameters at N 3 the equilibrium is unstable and a stable limit cycle appears from a Hopf bifurcation; for parameters at N 4 the equilibrium is an weak unstable focus Hopf point with l 1 positive and there is a stable limit cycle; for parameters at N 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 N 6 the equilibrium is asymptotically stable and the two cycles collide, giving rise to a nondegenerate fold bifurcation of the cycles.
Then, using the notation of Section 3 see expression 3.3 the multilinear symmetric functions corresponding to f can be written as

4.15
The eigenvalues of A are and from 3.5 one has q iω 0 ν 2μν iω 0 , iω 0 , β 2c , 4.17 where Using these calculations we prove the next three theorems, from which statement b of introduction follows. where
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 l 1 and the expression 3 2αμ − 1 ν 3 are negative.In Figure 2 b the surface {l 1 0} and the regions S and U are illustrated.These regions are defined implicitly as {l 1 < 0} and {l 1 > 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 Therefore, the transversality condition holds at the Hopf point.The theorem is proved.Phase portraits of system 1.1 for the flow restricted to the center manifold and its continuation.For parameters at M 1 the equilibrium is unstable; for parameters at M 2 the equilibrium is an weak unstable focus Hopf point with l 1 positive ; for parameters at M 3 the equilibrium is stable and an unstable limit cycle appears from a Hopf bifurcation; for parameters at M 4 the equilibrium is an weak stable focus Hopf point with l 1 negative and there is an unstable limit cycle; for parameters at M 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 M 6 the equilibrium is unstable and the two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles.
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 Figure 2 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 l 2 is positive.The graph of h μ, ν, α 0 on the surface {l 1 0}, which determines the curve C, and the signs of the first and second Lyapunov coefficients are also illustrated.As follows, l 2 < 0 on the open region on the surface {l 1 0}, denoted by C 1 .On this region a typical reference point R is depicted.Also l 2 > 0 on the open region on the surface {l 1 0}, denoted by C 2 .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.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 β β 2c , the solid curves where l 1 0, the regions S and U where l 1 < 0 and l 1 > 0, respectively, the dashed curves where l 2 0 and the arcs C 1 and C 2 where l 2 < 0 and l 2 > 0, respectively.The bifurcation diagrams for typical points R 1 and R 2 are depicted in Figure 3      coefficient at the point Q 1 , that is, for the case α 2 or equivalently for R R p .For the case α 3 the calculations are very similar.Theorem 4.6.Consider α 2 in system 1.1 .Then there is a unique point μ, ν, β 2c where the surfaces {l 1 0} and {l 2 0} intersect transversally on the Hopf hypersurface H 2 .For the parameter values at the point Q 1 , system 1.1 has a transversal Hopf point of codimension 3 at E which is unstable since l 3 Q 1 > 0. The bifurcation diagram of system 1.1 at the point Q 1 is illustrated in Figure 7.In Figure 8 the bifurcation diagram of system 1.1 at a typical point R 1 of Figure 7 is drawn.
, and V 9 in the parameter space near the point R 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 V i , i 1, . . ., 9 of Figure 8 a .

Computer Assisted Proof
The point Q 1 is the intersection of the surfaces {l

4.28
The previous calculations have also been corroborated with 20 decimals round-off precision performed using the software MATHEMATICA 5

4.29
The transversality condition at Q 1 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 8 > 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 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 l 2 < 0, l 1 > 0, and l 0 > 0 with |l 2 | l 1 l 0 > 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 V 1 the equilibrium is unstable; for parameters at V 2 the equilibrium is an weak unstable focus Hopf point with positive l 1 ; for parameters at V 3 the equilibrium is stable and an unstable limit cycle appears from the Hopf bifurcation; for parameters at V 4 the equilibrium is an weak stable focus Hopf point with negative l 1 and there is an unstable limit cycle; for parameters at V 5 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 V 6 the equilibrium is unstable and two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles; for parameters at V 7 the equilibrium is stable and two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles; for parameters at V 8 the equilibrium is stable and two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles; for parameters at V 9 the equilibrium is stable and there are three limit cycles encircling the equilibrium.

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 E 0 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 .

Bifurcations at E 0
For μ < 0, system 1.1 has the origin E 0 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 .
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 E 0 , which spiral to the equilibria E ± for The point E 0 is an unstable focus and a limit cycle arises around it.Eigenvalues at E 0 are λ 1 −2.774, λ 2,3 0.012 ± 0.728i.
β > β 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 .
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 E 0 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 .

Bifurcation at the Critical Point Q 1 in the Parameter Space
As stated in Theorem 4.6, the most complex Hopf bifurcation at the points E ± occurs at the point Q 1 μ, ν, β 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 .036345241, 4.023353871 .In Figure 16 we show an approximation of the unstable manifolds of the equilibrium E 0 .These manifolds spiral to the equilibria E ± if β > β Q 1 4.023353871 and this spiralling behavior become slower as β tends to β Q 1 β 2c from above for μ μ Q 1 and ν ν Q 1 fixed .
For μ, α, β at the point Q 1 see Theorem 4.6 the system presents a highly degenerate behavior.In fact, due to the vanishing of the first and second Lyapunov coefficients, l 1 and l 2 , 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 l 3 Q 1 > 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 Q 1 .
For β < β Q 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.

Concluding Remarks
In this paper we analyze the Lyapunov stability of the equilibria E 0 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 E 0 resp., E ± is pushed forward to the calculation of the second resp., second and third .The three equilibria are unstable and there exists a "big" stable periodic orbit encircling them.Eigenvalues at E 0 are λ 1 0.368, λ 2 −1.969, λ 3 −0.337.Eigenvalues at E ± are λ 1 −2.144, λ 2,3 0.010 ± 0.478i.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 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 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.

Figure 1 :
Figure 1: ADVP circuit with parallel resistor R p .

Theorem 4 . 1 .
Consider the four-parameter family of differential equations 1.1 .The first Lyapunov coefficient associated to the equilibrium E 0 is given by

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

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 l 2 associated to the equilibrium E is given by

Figure 2 (
b)), then system 1.1 has a transversal Hopf point of codimension 2 at E .More specifically, if μ, ν, α, β 2c ∈ C 1 then the Hopf point at E is asymptotically stable and the bifurcation diagram is illustrated in Figure3for a typical point R.If μ, ν, α, β 2c ∈ C 2 then the Hopf point at E is unstable and the bifurcation diagram is drawn in Figure4for a typical point T .

Figure 4 :
Figure4: a Typical points M 1 , M 2 , M 3 , M 4 , M 5 , and M 6 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 M 1 the equilibrium is unstable; for parameters at M 2 the equilibrium is an weak unstable focus Hopf point with l 1 positive ; for parameters at M 3 the equilibrium is stable and an unstable limit cycle appears from a Hopf bifurcation; for parameters at M 4 the equilibrium is an weak stable focus Hopf point with l 1 negative and there is an unstable limit cycle; for parameters at M 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 M 6 the equilibrium is unstable and the two cycles collide corresponding to a nondegenerate fold bifurcation of the cycles.

Remark 4 . 5 .
There are numerical evidences that the sign of the third Lyapunov coefficient is always positive along the curve C given by the intersection of the surfaces {l 1 0} and {l 2 0} see Figure2 b.For α 1 see the article in 5 .For α 2 see Theorem 4.6 in what follows.Figure5shows the behavior of the third Lyapunov coefficient as a function of α, for μ, ν, β on the curve C.So the bifurcation diagrams depicted in Figures7 and 8are essentially the same for all values of α ≥ 1.

Figure 5 :
Figure 5: Graph of l 3 with respect to α.