# On the Computation of Degenerate Hopf Bifurcations forn-Dimensional Multiparameter Vector Fields

Authors: Michail P. Markakis; Panagiotis S. Douris
Journal: International Journal of Mathematics and Mathematical Sciences
 (2016)

---

## Abstract

The restriction of ann-dimensional nonlinear parametric system on the center manifold is treated via a new proper symbolic form and analytical expressions of the involved quantities are obtained as functions of the parameters by lengthy algebraic manipulations combined with computer assisted calculations. Normal forms regarding degenerate Hopf bifurcations up to codimension 3, as well as the corresponding Lyapunov coefficients and bifurcation portraits, can be easily computed for any system under consideration.


---

## Body

## 1. Introduction

In the international literature, we encounter a large number ofn-dimensional nonlinear dynamical systems (n≥3) in different areas of applied sciences, which are studied with respect to stability and local bifurcations generated by the variation of the involved parameters. We refer, for example, to recent relevant works in modern scientific disciplines like biosciences [1–4], energy systems [5, 6], economics [7, 8], and so forth. We also refer to classic dynamical systems like the Lorenz [9], Lü [10, 11], and Chen [12] systems or to more general modified forms like [13].Moreover, regarding the Hopf bifurcation, the combination of the center manifold theory with the Poincaré normal forms for planar systems, applied ton-dimensional ones, is analytically presented in the classic book of Kuznetsov [14] and the relevant references therein. The method is based on the state space decomposition in the critical and noncritical subspace in addition to the necessary normalization with respect to the involved eigenvectors. Then substitution of the reduced critical equations in the invariance relation concerning the center manifold (Wc) ([14], Sections 5.4.1, 8.7.1, and 8.7.3) leads to the computation of the critical polynomial coefficients of Wc and therefore to the evaluation of the coefficients of the equation restricted to the center manifold. Finally, according to the obtained two-dimensional normal forms ([14], Sections 3.5 and 8.3), formulae for the resonant odd terms and the corresponding critical Lyapunov coefficients are obtained up to codimension 2 (Bautin bifurcation), and the corresponding bifurcation diagram is also shown. In other works, like that of Sotomayor et al. [15], based on the same procedure, critical Lyapunov coefficients are extracted for degenerate Hopf bifurcations up to codimension 4. The same authors deal with a three-parameter system, obtaining results with respect to degenerate cases up to codimension 3 [16].Regarding the present work, on the one hand, the paper stands for a review of the analytical procedure concerning the restriction of ann-dimensional nonlinear three-parameter system on the center manifold, as well as the normal forms of a planar system with respect to degenerate Hopf bifurcations of higher codimension. More precisely, as regards the normal transformations of an equation in the zz- plane, we obtain the locally orbitally equivalent systems and the corresponding expressions of Lyapunov coefficients up to codimension 3. Moreover the bifurcation portraits of the codimension 2 and 3 cases are presented in detail in Section 3 (Section 3.2).On the other hand, a new appropriate symbolic representation is adopted, by which we easily attain specification of the necessary (according to the codimension under consideration) order of the center manifold that should be kept in each term involved in the analytic calculations referring to the derivation of the restricted equation and the treatment of the invariance relation, as well. Hence considering the required extensive algebraic manipulations and taking into account the part assigned to the computer (carried out by use of the symbolic algebra software Mathematica 7), the procedure used herein (Sections2, 2.2, and 2.3) is proved excellent at saving computational time and memory space.Furthermore we note that the state space decomposition is applied here on the equilibrium path and hence the reduced equations and all the formulae obtained are expressed in terms of the parameters of the system. In particular the center manifold coefficients are derived as parameter functions explicitly and implicitly via the eigen-quantities, as well as the lower order coefficients ofWc. The obtained lengthy (in spite of the implicit structure) expressions up to the sixth order are stored in Mathematica files, ready to be computed throughout the region of the parameter space of interest. Therefore the coefficients of the restricted equation expressed in terms of the eigenvectors and the center manifold coefficients and also theresonant terms of the considered normal forms (being functions of the planar (restricted) coefficients), as well as the corresponding Lyapunov coefficients, are numerically evaluated fast, not only at the critical parameter values, but also throughout the whole parameter space. Thus the respective bifurcation portraits can be constructed. Moreover, the parameter-dependent formulae allow the computer assisted calculation of the derivatives of the various coefficients involved in the analysis and so based on the same implicit structure, thetransversality condition corresponding to the codimension of the considered bifurcation can be easily verified at the critical equilibrium. The basic steps of the followed procedure are displayed in Section 2.3.We note that the analysis concerns general multiparameter systems (we refer to the three-parameter case herein), in order to take into account the combined effects generated by the variation of as many as possible significant quantities of the physical background of the problem. We should also note that the form and the structure of the representation adopted, as well as the obtained analytic formulae, allow the generalization and application of this procedure in any multidimensional and multiparameter system, with minor modifications. Finally we would like to highlight that, as mentioned before, it is the structure and the steps of the analytical procedure that are shown here, while the lengthy formulae obtained, not listed here, are included in the respective electronic files, produced via algorithms based on the symbolic algebraic form of the analysis adopted herein.


## 2. Analysis and Formulae for the Center Manifold Theory

### 2.1. Reduction to an(n+2)-Dimensional Coordinate Space

Consider a smooth continuous-time three-parameter system with smooth dependence on the parameters:(1)dxdt=fx;awith x∈Rn, a=(a1,a2,a3)T∈R3, and f:Rn+3→Rn with f∈C∞. If x0(a) is the equilibrium path of (1) and J0(a)=Dxf(x0(a);a), the Jacobian matrix evaluated at this equilibrium path, then, for small perturbations ξ=x-x0(a), expansion of (1) yields(2)dξdt=J0aξ+Fξ;a.The smooth vector function F:Rn+3→Rn represents the nonlinear terms of the right-hand side of (2); that is, F=Ox2 and(3)Fξ;a=12F2ξ,ξ;a+16F3ξ,ξ,ξ;a+⋯,where(4)Fmu1,…,um;a=∑j1,…,jm=1n∂mfx;a∂xj1⋯∂xjmx=x0uj11⋯ujmm,m=2,3,…with u1,…,um∈Cn. It is well known that the Hopf bifurcation theorem (see [17], Section 11.2) implies that if, at a critical triplet a0=a10,a20,a30, J0 has a pair of pure imaginary eigenvalues λ1,2=±iω0, ω0>0, while the real part of the remaining eigenvalues is negative, then an oscillatory type of instability occurs close to a0, leading to a family of limit cycles. Moreover we conclude that, in a small region around a0, J0 has a pair of complex conjugate eigenvalues λa and λ-a with(5)λa=μa+iωa,μa0=0,ωa0=ω0>0.We also note that the classic Hopf theory additionally demands that, at a0, λ(a) cross the imaginary axis with nonzero speed (transversality condition).Let us now consider the complex eigenvectorsqa∈Cn and pa∈Cn of J0 and J0T, respectively, having the properties(6)J0q=λq,J0Tp=λ-pand satisfying the normalization condition(7)p,q=∑i=1np-iqi=1.By the above normalization it can also be easily proved that(8)p,q-=0.Let T(α) denote the two-dimensional parameter-dependent real eigenspace (corresponding to λ,λ-), spanned by Re⁡[qα],Im⁡[qα], which becomes critical at a0(Tc=T(α0)), and Tsu(α) the (n-2)-dimensional (n>2) real eigenspace corresponding to all eigenvalues of J0 other than λ,λ-. Then, by use of the following Lemma ([14], Lemma 5.4) we get the following.Lemma 1.
Considery∈Tsu if and only if 〈p,y〉=0.Considery∈Rn; we introduce a complex variable z∈C1 and decompose any x∈Rn as(9)x-x0=ξ=zq+z-q-+y,where zq+z-q-∈T and y∈Tsu. Taking into account (7) and (8), the last equation yields(10)z=p,ξ,y=ξ-p,ξq-p-,ξq-.Thus taking into account (2), as well as (6) and (9), differentiation of (10) yields the reduced form of system (2) in the (n+2)-dimensional coordinate space (z,y):(11)z˙=λaz+gz,z-,y;a,(12)y˙=J0ay+Fzq+z-q-+y;a-gz,z-,y;aq-g-z,z-,y;aq-,where y˙=dy/dt and(13)gz,z-,y;a=p,Fzq+z-q-+y;a.Systems (2) and (11)-(12) are in fact dimensionally equivalent, due to the two real constraints imposed on y by Lemma 1.


### 2.2. A New Symbolic Representation of System (11)-(12)

The complex functiong(z,z-), given by (13), has the following form:(14)g=12g002+16g003+⋯+g101+12g102+16g103+⋯z+g011+12g012+16g013+⋯z-+g2+g3+⋯,where(15)gmz,z-,y=∑k+l=m1k!l!gkl+gkl1+12gkl2+16gkl3+⋯zkz-l,gijm+y=p,Fi+j+m+q,…,q︸i,q-,…,q-︸j,y,…,y︸m+,gkl=gkl0with gij(m+)∈C, m=2,3,…, i+j=0,1,2,…, m+=0,1,2,….Thus by using (14), (11) and (12) become(16)z˙=12g002+16g003+⋯+λ+g101+12g102+16g103+⋯z+g011+12g012+16g013+⋯z-+g2+g3+⋯,(17)y˙=J0ay+12r002+16r003+⋯-12g002+16g003+⋯q-12g-002+16g-003+⋯q-+r101+12r102+16r103+⋯-g101+12g102+16g103+⋯q-g-011+12g-012+16g-013+⋯q-z+r011+12r012+16r013+⋯-g011+12g012+16g013+⋯q-g-101+12g-102+16g-103+⋯q-z-+r2+r3+⋯,where(18)rmz,z-,y=∑k+l=m1k!l!rkl+rkl1+12rkl2+16rkl3+⋯-gkl+gkl1+12gkl2+16gkl3+⋯q-g-lk+g-lk1+12g-lk2+16g-lk3+⋯q-zkz-l,rijm+y=Fi+j+m+q,…,q︸i,q-,…,q-︸j,y,…,y︸m+,rkl=rkl0with gij(m+) as in (15), rij(m+)∈Cn, m=2,3,…, i+j=0,1,2,…, m+=0,1,2,….


### 2.3. Computation of the Center Manifold

Regarding now the center manifoldWac, we adopt the representation(19)y=Vz,z-=w2z,z-+w3z,z-+⋯,wmz,z-=∑k+l=m1k!l!wklzkz-l,m=2,3,…with wkl∈Cn, k+l=2,3,…, and 〈p,wkl〉=0. Since we want to investigate the system under consideration up to codimension 3 bifurcation of the equilibrium, by substituting in (16) for y the center manifold expression (19) and keeping terms up to seventh order with respect to z, we arrive at(20)z˙=12Q2,0,0,2,0,0,0,0+12Q2,0,0,0,2,0,0,0+Q2,0,0,1,1,0,0,0+Q2,0,0,1,0,1,0,0+Q2,0,0,1,0,0,1,0+Q2,0,0,0,1,1,0,0+16Q3,0,0,3,0,0,0,0+12Q3,0,0,2,1,0,0,0+λz+∑k+l=1Q2,k,l,1,0,0,0,0+Q2,k,l,0,1,0,0,0+Q2,k,l,0,0,1,0,0+Q2,k,l,0,0,0,1,0+Q2,k,l,0,0,0,0,1+12Q3,k,l,2,0,0,0,0+12Q3,k,l,0,2,0,0,0+Q3,k,l,1,1,0,0,0+Q3,k,l,1,0,1,0,0+16Q4,k,l,3,0,0,0,0zkz-l+∑k+l=21k!l!gkl+Q3,k,l,1,0,0,0,0+Q3,k,l,0,1,0,0,0+Q3,k,l,0,0,1,0,0+Q3,k,l,0,0,0,1,0+12Q4,k,l,2,0,0,0,0+Q4,k,l,1,1,0,0,0zkz-l+∑k+l=31k!l!gkl+Q4,k,l,1,0,0,0,0+Q4,k,l,0,1,0,0,0+Q4,k,l,0,0,1,0,0+12Q5,k,l,2,0,0,0,0zkz-l+∑k+l=41k!l!gkl+Q5,k,l,1,0,0,0,0+Q5,k,l,0,1,0,0,0zkz-l+∑k+l=51k!l!gkl+Q6,k,l,1,0,0,0,0zkz-l+1720g60z6+1120g51z5z-+148g42z4z-2+136g33z3z-3+148g24z2z-4+1120g15zz-5+1720g06z-6+15040g70z7+1720g61z6z-+1240g52z5z-2+1144g43z4z-3+1144g34z3z-4+1240g25z2z-5+1720g16zz-6+15040g07z-7+Oz8,where(21)Qm,i,j,m2+,m3+,m4+,m5+,m6+=p,Hm,i,j,m2+,m3+,m4+,m5+,m6+,Hm,i,j,m2+,m3+,m4+,m5+,m6+=Fmq,…,q︸i,q-,…,q-︸j,w2,…,w2︸m2+,w3,…,w3︸m3+,w4,…,w4︸m4+,w5,…,w5︸m5+,w6,…,w6︸m6+with m=i+j+m2++m3++m4++m5++m6+. By means of (20) we conclude that we need to evaluate the center manifold coefficients up to the sixth order (wkl, k+l=2,…,6). Thus by writing the invariance relation for Wac,(22)y˙=Vzz˙+Vz-z-˙=W1+W2+W3+W4+W5+⋯z˙+W~1+W~2+W~3+W~4+W~5+⋯z-˙,where(23)Wmz,z-=∂wm+1z,z-∂z,W~mz,z-=∂wm+1z,z-∂z-,m=1,2,…,we substitute (16) and (17) in (22), where the first equation of (19) is further substituted for y in gijm+y and rijm+y, i+j=0,1,2,…, m+=1,2,… (see (15) and (18)), and keep terms up to the sixth order with respect to z. Then, by considering the resulting equation, by means of tedious algebra combined with calculations based on symbolic computation software (here Mathematica 7), we proceed as figured out in the following steps:(1)
By using notation (21), that is, H, Q and Q- for rijm+w(2)+w(3)+⋯, gijm+w(2)+w(3)+⋯, and g-ijm+w(2)+w(3)+⋯, respectively, we equate the terms of the same order, from the second up to the sixth and obtain five equations, namely, Ei, i=2,…,6.(2)
In the obtained equations(Ei), by substituting the second equation of (19) for w(m), m=2,…,6, as well as the associated expressions giving their derivatives with respect to z and z-, for W(m-1) and W~(m-1), m=2,…,6, respectively, we conclude with new explicit forms of the above equations, namely, (Ei′), i=2,…,6.(3)
In each equation(Ei′), i=2,…,6, by considering the coefficients of the zkz-l terms, k+l=i, and taking into account the following relations,(24)J0-nλ+mλ-I-1q=-1n-1λ+mλ-q,J0-nλ+mλ-I-1q-=-1nλ+m-1λ-q-,based on the first equation of (6), with n+m=2,…,6, we solve for wkl:(25)wkl=wkla,μa,ωa,qa,q-a,pa,p-a,wij,w-ij,i+j=2,…,k+l-1.Considerk+l=2,…,6. In particular in the critical case μ=0, ω=ω0, for k-l=±1, since the coefficient matrix of the corresponding system for wkl becomes J0∓iω0I with determinant equal to zero, the solution is obtained in combination with the additional condition 〈p,wkl〉=0 (see Section 2.4).(4)
By substituting the above derived expressions ofwkl, k+l=2,…,6, in (20) and keeping terms up to seventh order with respect to z, we obtain the coefficients Gkl, k+l=2,…,7, of the two-dimensional equation in the zz- plane, representing the restriction of the system on the center manifold; namely,(26)z˙=λaz+G2+⋯+G7+Oz8,Gmz,z-=∑k+l=m1k!l!Gklzkz-l,m=2,…,7,Gkl=Gklqa,q-a,p-a,k+l=2,Gkl=Gklqa,q-a,p-a,wij,k+l≥3,i+j=2,…,k+l-1.(5)
Finally, by considering the analytic expressions (39) of theresonant terms ci corresponding to the respective normal forms (i=1,2, or 3, depending on the considered bifurcation codimension; see Section 3.1 below), we substitute gi+1,i, gst, and g-st(s+t=2,…,2i) with the coefficients Gi+1,i, Gst, and G-st, respectively, of the equation restricted to the center manifold. Then the necessary Lyapunov coefficients (provided in Section 3.1 below) can be computed numerically fast by any reliable symbolic mathematical computation software, not only at the critical parameter values, but also throughout the region of the parameter space of interest. Hence the bifurcation portraits of degenerate cases (given in Section 3.2) can be easily constructed. Moreover, the induced evaluation of the critical derivatives concludes with the verification of the correspondingtransversality conditions (see C1-2, C2-2, and (C3-2a, C3-2b, and C3-2c) in Sections 3.1.1, 3.1.2, and 3.1.3, resp.).


### 2.4. Derivation of the Critical Center Manifold Coefficientswkl, k-l=±1

Due to the form of the equations resulting in the vector coefficientswkl, namely,(27)J0-kλ+lλ-Iwkl=Akl,λ=μ+iω,or(28)J0-k+lμ+ik-lωIwkl=Akl,for k-l=±1, in the critical case μ=0, ω=ω0, the determinant becomes zero. Here by means of the constraint(29)p,wkl=0,at first it can be easily proved that 〈p,J0∓iωIwkl〉=0 and hence 〈p,Akl〉=0, which, as it can further be proved, is a sufficient and necessary condition for the solvability of the corresponding system; namely,(30)J0∓iω0Iwkl=Akl,k-l=±1.The expressions of the solutions for the critical coefficients wkl,k-l=±1, also satisfy constraint (29). Moreover as regards the expressions of the critical derivatives of wkl,k-l=±1, with respect to aι,ι=1,2,3, involved in thetransversality condition, we solve the system obtained by differentiation of (28), taking into account the respective derivative of (29).Thus aiming to form a compact notation introduced in symbolic calculations, based on (30) and (29), as well as their parametric derivatives, we write(31)X=x1A-Bx2A-x1KAD33+BD32-ΓD33AD33Z,Y=x2A-x1KD33-D32D33Z,Z=εAD33+p-1Bx2A-x1K-x1D33-p-2Ax2A-x1Kp-1BD32-ΓD33-p-2AD32+p-3AD33,where D32, D33 are minor determinants of the matrix:(32)
J
0
∓
i
ω
0
I
=
A
B
Γ
K
Λ
M
ε
1
ε
2
ε
3
,
k
-
l
=
±
1
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
X
,
Y
,
Z
x
1
,
x
2
,
x
3
ε
w
k
l
A
k
l
0
d
w
k
l
d
a
i
d
A
k
l
d
a
i
-
d
J
0
d
a
i
-
k
+
l
d
μ
d
a
i
+
i
k
-
l
d
ω
d
a
i
I
w
k
l
=
B
k
l
i
,
-
∑
n
=
1
3
d
p
-
n
d
a
i
w
k
l
nHere, we also provide the expressions of the vectors q and p:(33)q=L31,-L32,L33with L3i, i=1,2,3, minor determinants of J0-μ+iωI and(34)p=1L31Tq-1-L32Tq-2+L33Tq-3L31T,-L32T,L33Twith L3iT, i=1,2,3, minor determinants of J0T-μ-iωI.


## 2.1. Reduction to an(n+2)-Dimensional Coordinate Space

Consider a smooth continuous-time three-parameter system with smooth dependence on the parameters:(1)dxdt=fx;awith x∈Rn, a=(a1,a2,a3)T∈R3, and f:Rn+3→Rn with f∈C∞. If x0(a) is the equilibrium path of (1) and J0(a)=Dxf(x0(a);a), the Jacobian matrix evaluated at this equilibrium path, then, for small perturbations ξ=x-x0(a), expansion of (1) yields(2)dξdt=J0aξ+Fξ;a.The smooth vector function F:Rn+3→Rn represents the nonlinear terms of the right-hand side of (2); that is, F=Ox2 and(3)Fξ;a=12F2ξ,ξ;a+16F3ξ,ξ,ξ;a+⋯,where(4)Fmu1,…,um;a=∑j1,…,jm=1n∂mfx;a∂xj1⋯∂xjmx=x0uj11⋯ujmm,m=2,3,…with u1,…,um∈Cn. It is well known that the Hopf bifurcation theorem (see [17], Section 11.2) implies that if, at a critical triplet a0=a10,a20,a30, J0 has a pair of pure imaginary eigenvalues λ1,2=±iω0, ω0>0, while the real part of the remaining eigenvalues is negative, then an oscillatory type of instability occurs close to a0, leading to a family of limit cycles. Moreover we conclude that, in a small region around a0, J0 has a pair of complex conjugate eigenvalues λa and λ-a with(5)λa=μa+iωa,μa0=0,ωa0=ω0>0.We also note that the classic Hopf theory additionally demands that, at a0, λ(a) cross the imaginary axis with nonzero speed (transversality condition).Let us now consider the complex eigenvectorsqa∈Cn and pa∈Cn of J0 and J0T, respectively, having the properties(6)J0q=λq,J0Tp=λ-pand satisfying the normalization condition(7)p,q=∑i=1np-iqi=1.By the above normalization it can also be easily proved that(8)p,q-=0.Let T(α) denote the two-dimensional parameter-dependent real eigenspace (corresponding to λ,λ-), spanned by Re⁡[qα],Im⁡[qα], which becomes critical at a0(Tc=T(α0)), and Tsu(α) the (n-2)-dimensional (n>2) real eigenspace corresponding to all eigenvalues of J0 other than λ,λ-. Then, by use of the following Lemma ([14], Lemma 5.4) we get the following.Lemma 1.
Considery∈Tsu if and only if 〈p,y〉=0.Considery∈Rn; we introduce a complex variable z∈C1 and decompose any x∈Rn as(9)x-x0=ξ=zq+z-q-+y,where zq+z-q-∈T and y∈Tsu. Taking into account (7) and (8), the last equation yields(10)z=p,ξ,y=ξ-p,ξq-p-,ξq-.Thus taking into account (2), as well as (6) and (9), differentiation of (10) yields the reduced form of system (2) in the (n+2)-dimensional coordinate space (z,y):(11)z˙=λaz+gz,z-,y;a,(12)y˙=J0ay+Fzq+z-q-+y;a-gz,z-,y;aq-g-z,z-,y;aq-,where y˙=dy/dt and(13)gz,z-,y;a=p,Fzq+z-q-+y;a.Systems (2) and (11)-(12) are in fact dimensionally equivalent, due to the two real constraints imposed on y by Lemma 1.


## 2.2. A New Symbolic Representation of System (11)-(12)

The complex functiong(z,z-), given by (13), has the following form:(14)g=12g002+16g003+⋯+g101+12g102+16g103+⋯z+g011+12g012+16g013+⋯z-+g2+g3+⋯,where(15)gmz,z-,y=∑k+l=m1k!l!gkl+gkl1+12gkl2+16gkl3+⋯zkz-l,gijm+y=p,Fi+j+m+q,…,q︸i,q-,…,q-︸j,y,…,y︸m+,gkl=gkl0with gij(m+)∈C, m=2,3,…, i+j=0,1,2,…, m+=0,1,2,….Thus by using (14), (11) and (12) become(16)z˙=12g002+16g003+⋯+λ+g101+12g102+16g103+⋯z+g011+12g012+16g013+⋯z-+g2+g3+⋯,(17)y˙=J0ay+12r002+16r003+⋯-12g002+16g003+⋯q-12g-002+16g-003+⋯q-+r101+12r102+16r103+⋯-g101+12g102+16g103+⋯q-g-011+12g-012+16g-013+⋯q-z+r011+12r012+16r013+⋯-g011+12g012+16g013+⋯q-g-101+12g-102+16g-103+⋯q-z-+r2+r3+⋯,where(18)rmz,z-,y=∑k+l=m1k!l!rkl+rkl1+12rkl2+16rkl3+⋯-gkl+gkl1+12gkl2+16gkl3+⋯q-g-lk+g-lk1+12g-lk2+16g-lk3+⋯q-zkz-l,rijm+y=Fi+j+m+q,…,q︸i,q-,…,q-︸j,y,…,y︸m+,rkl=rkl0with gij(m+) as in (15), rij(m+)∈Cn, m=2,3,…, i+j=0,1,2,…, m+=0,1,2,….


## 2.3. Computation of the Center Manifold

Regarding now the center manifoldWac, we adopt the representation(19)y=Vz,z-=w2z,z-+w3z,z-+⋯,wmz,z-=∑k+l=m1k!l!wklzkz-l,m=2,3,…with wkl∈Cn, k+l=2,3,…, and 〈p,wkl〉=0. Since we want to investigate the system under consideration up to codimension 3 bifurcation of the equilibrium, by substituting in (16) for y the center manifold expression (19) and keeping terms up to seventh order with respect to z, we arrive at(20)z˙=12Q2,0,0,2,0,0,0,0+12Q2,0,0,0,2,0,0,0+Q2,0,0,1,1,0,0,0+Q2,0,0,1,0,1,0,0+Q2,0,0,1,0,0,1,0+Q2,0,0,0,1,1,0,0+16Q3,0,0,3,0,0,0,0+12Q3,0,0,2,1,0,0,0+λz+∑k+l=1Q2,k,l,1,0,0,0,0+Q2,k,l,0,1,0,0,0+Q2,k,l,0,0,1,0,0+Q2,k,l,0,0,0,1,0+Q2,k,l,0,0,0,0,1+12Q3,k,l,2,0,0,0,0+12Q3,k,l,0,2,0,0,0+Q3,k,l,1,1,0,0,0+Q3,k,l,1,0,1,0,0+16Q4,k,l,3,0,0,0,0zkz-l+∑k+l=21k!l!gkl+Q3,k,l,1,0,0,0,0+Q3,k,l,0,1,0,0,0+Q3,k,l,0,0,1,0,0+Q3,k,l,0,0,0,1,0+12Q4,k,l,2,0,0,0,0+Q4,k,l,1,1,0,0,0zkz-l+∑k+l=31k!l!gkl+Q4,k,l,1,0,0,0,0+Q4,k,l,0,1,0,0,0+Q4,k,l,0,0,1,0,0+12Q5,k,l,2,0,0,0,0zkz-l+∑k+l=41k!l!gkl+Q5,k,l,1,0,0,0,0+Q5,k,l,0,1,0,0,0zkz-l+∑k+l=51k!l!gkl+Q6,k,l,1,0,0,0,0zkz-l+1720g60z6+1120g51z5z-+148g42z4z-2+136g33z3z-3+148g24z2z-4+1120g15zz-5+1720g06z-6+15040g70z7+1720g61z6z-+1240g52z5z-2+1144g43z4z-3+1144g34z3z-4+1240g25z2z-5+1720g16zz-6+15040g07z-7+Oz8,where(21)Qm,i,j,m2+,m3+,m4+,m5+,m6+=p,Hm,i,j,m2+,m3+,m4+,m5+,m6+,Hm,i,j,m2+,m3+,m4+,m5+,m6+=Fmq,…,q︸i,q-,…,q-︸j,w2,…,w2︸m2+,w3,…,w3︸m3+,w4,…,w4︸m4+,w5,…,w5︸m5+,w6,…,w6︸m6+with m=i+j+m2++m3++m4++m5++m6+. By means of (20) we conclude that we need to evaluate the center manifold coefficients up to the sixth order (wkl, k+l=2,…,6). Thus by writing the invariance relation for Wac,(22)y˙=Vzz˙+Vz-z-˙=W1+W2+W3+W4+W5+⋯z˙+W~1+W~2+W~3+W~4+W~5+⋯z-˙,where(23)Wmz,z-=∂wm+1z,z-∂z,W~mz,z-=∂wm+1z,z-∂z-,m=1,2,…,we substitute (16) and (17) in (22), where the first equation of (19) is further substituted for y in gijm+y and rijm+y, i+j=0,1,2,…, m+=1,2,… (see (15) and (18)), and keep terms up to the sixth order with respect to z. Then, by considering the resulting equation, by means of tedious algebra combined with calculations based on symbolic computation software (here Mathematica 7), we proceed as figured out in the following steps:(1)
By using notation (21), that is, H, Q and Q- for rijm+w(2)+w(3)+⋯, gijm+w(2)+w(3)+⋯, and g-ijm+w(2)+w(3)+⋯, respectively, we equate the terms of the same order, from the second up to the sixth and obtain five equations, namely, Ei, i=2,…,6.(2)
In the obtained equations(Ei), by substituting the second equation of (19) for w(m), m=2,…,6, as well as the associated expressions giving their derivatives with respect to z and z-, for W(m-1) and W~(m-1), m=2,…,6, respectively, we conclude with new explicit forms of the above equations, namely, (Ei′), i=2,…,6.(3)
In each equation(Ei′), i=2,…,6, by considering the coefficients of the zkz-l terms, k+l=i, and taking into account the following relations,(24)J0-nλ+mλ-I-1q=-1n-1λ+mλ-q,J0-nλ+mλ-I-1q-=-1nλ+m-1λ-q-,based on the first equation of (6), with n+m=2,…,6, we solve for wkl:(25)wkl=wkla,μa,ωa,qa,q-a,pa,p-a,wij,w-ij,i+j=2,…,k+l-1.Considerk+l=2,…,6. In particular in the critical case μ=0, ω=ω0, for k-l=±1, since the coefficient matrix of the corresponding system for wkl becomes J0∓iω0I with determinant equal to zero, the solution is obtained in combination with the additional condition 〈p,wkl〉=0 (see Section 2.4).(4)
By substituting the above derived expressions ofwkl, k+l=2,…,6, in (20) and keeping terms up to seventh order with respect to z, we obtain the coefficients Gkl, k+l=2,…,7, of the two-dimensional equation in the zz- plane, representing the restriction of the system on the center manifold; namely,(26)z˙=λaz+G2+⋯+G7+Oz8,Gmz,z-=∑k+l=m1k!l!Gklzkz-l,m=2,…,7,Gkl=Gklqa,q-a,p-a,k+l=2,Gkl=Gklqa,q-a,p-a,wij,k+l≥3,i+j=2,…,k+l-1.(5)
Finally, by considering the analytic expressions (39) of theresonant terms ci corresponding to the respective normal forms (i=1,2, or 3, depending on the considered bifurcation codimension; see Section 3.1 below), we substitute gi+1,i, gst, and g-st(s+t=2,…,2i) with the coefficients Gi+1,i, Gst, and G-st, respectively, of the equation restricted to the center manifold. Then the necessary Lyapunov coefficients (provided in Section 3.1 below) can be computed numerically fast by any reliable symbolic mathematical computation software, not only at the critical parameter values, but also throughout the region of the parameter space of interest. Hence the bifurcation portraits of degenerate cases (given in Section 3.2) can be easily constructed. Moreover, the induced evaluation of the critical derivatives concludes with the verification of the correspondingtransversality conditions (see C1-2, C2-2, and (C3-2a, C3-2b, and C3-2c) in Sections 3.1.1, 3.1.2, and 3.1.3, resp.).


## 2.4. Derivation of the Critical Center Manifold Coefficientswkl, k-l=±1

Due to the form of the equations resulting in the vector coefficientswkl, namely,(27)J0-kλ+lλ-Iwkl=Akl,λ=μ+iω,or(28)J0-k+lμ+ik-lωIwkl=Akl,for k-l=±1, in the critical case μ=0, ω=ω0, the determinant becomes zero. Here by means of the constraint(29)p,wkl=0,at first it can be easily proved that 〈p,J0∓iωIwkl〉=0 and hence 〈p,Akl〉=0, which, as it can further be proved, is a sufficient and necessary condition for the solvability of the corresponding system; namely,(30)J0∓iω0Iwkl=Akl,k-l=±1.The expressions of the solutions for the critical coefficients wkl,k-l=±1, also satisfy constraint (29). Moreover as regards the expressions of the critical derivatives of wkl,k-l=±1, with respect to aι,ι=1,2,3, involved in thetransversality condition, we solve the system obtained by differentiation of (28), taking into account the respective derivative of (29).Thus aiming to form a compact notation introduced in symbolic calculations, based on (30) and (29), as well as their parametric derivatives, we write(31)X=x1A-Bx2A-x1KAD33+BD32-ΓD33AD33Z,Y=x2A-x1KD33-D32D33Z,Z=εAD33+p-1Bx2A-x1K-x1D33-p-2Ax2A-x1Kp-1BD32-ΓD33-p-2AD32+p-3AD33,where D32, D33 are minor determinants of the matrix:(32)
J
0
∓
i
ω
0
I
=
A
B
Γ
K
Λ
M
ε
1
ε
2
ε
3
,
k
-
l
=
±
1
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
X
,
Y
,
Z
x
1
,
x
2
,
x
3
ε
w
k
l
A
k
l
0
d
w
k
l
d
a
i
d
A
k
l
d
a
i
-
d
J
0
d
a
i
-
k
+
l
d
μ
d
a
i
+
i
k
-
l
d
ω
d
a
i
I
w
k
l
=
B
k
l
i
,
-
∑
n
=
1
3
d
p
-
n
d
a
i
w
k
l
nHere, we also provide the expressions of the vectors q and p:(33)q=L31,-L32,L33with L3i, i=1,2,3, minor determinants of J0-μ+iωI and(34)p=1L31Tq-1-L32Tq-2+L33Tq-3L31T,-L32T,L33Twith L3iT, i=1,2,3, minor determinants of J0T-μ-iωI.


## 3. Normal Forms and Degenerate Bifurcations

### 3.1. Planar Case(n=2): Poincaré Normal Forms, Codimension nr

In this case, the termsgijm+ with i+j=0,1,2,…, m+=1,2,… in the first equation of (15) and in (16) vanish. Moreover, (17) does not exist. Now, in order to obtain the desired normal form, we first introduce the transformation(35)z=w+h2w,w-+⋯+h2nr+1w,w-,nr∈N,hmw,w-=∑k+l=m1k!l!hklwkw-l,m=2,…,2nr+1,the inversion of which is given by (see [14])(36)w=z-h2z,z--⋯-h2nr+1z,z-+Oz2nr+2.Then by differentiating (36) with respect to t and substituting the right-hand side of (11) for z˙ (with g=g2+g3+⋯), as well as the corresponding conjugate equation for z-˙, by keeping terms up to 2nr+1 order with respect to z, we obtain(37)w˙=λz+g2+⋯+g2nr+1-∑m=12nr-1λz+g2+⋯+g2nr+1-mHm-λ-z-+g-2+⋯+g-2nr+1-mH~m-λzH2nr-λ-z-H~2nr+Oz2nr+2,Hmz,z-=∂hm+1z,z-∂z,H~mz,z-=∂hm+1z,z-∂z-,m=1,…,2nr.Finally substitution of transformation (35) for z (as well as the respective conjugate relation for z-) in (37) yields(38)w˙=λaw+K2+⋯+K2nr+1+Ow2nr+2,Kmw,w-;a=∑k+l=mKklλ,λ-,hkl,gkl,hij,h-ij,gij,g-ijwkw-l,m=2,…,2nr+1,i+j=2,…,m-1,where terms up to 2nr+1 order with respect to w have been evaluated. Now, by setting Kkl, k+l=2,…,2nr+1, equal to zero we define the coefficients hkl of the transformation (35), while the remainingresonant terms of the odd order, ci=Ki+1,i, i=1,…,nr, constitute the desirable normal form of (11). Furthermore, by substituting the obtained expressions of hkl,h-kl involved in ci (k+l=2,…,2i), then by means of computer calculations we arrive at the analytic formulae of the normal coefficients:(39)ci=ciλ,λ-,gi+1,i,gst,g-st,s+t=2,…,2i,i=1,…,nr,which are lengthy for i≥2.


#### 3.1.1. Normal Form of the Third Order(nr=1)

In this case, the above procedure results in the third-ordernormal form:(40)w˙=λaw+c1aww2+Ow4,where c1(a)w2w- represents the cubicresonant term, due to the term (λ+λ-)h21 included in K21 (we set h21=0). Then by means of the linear time scaling(41)θ=ωatand the nonlinear time reparametrization(42)dθ=1-Imc1aωaw2dτ, (40) becomes(43)dwdτ=β+iw+l1aww2+Ow4,βa=μaωa,l1a=Rec1aωa-βaImc1aωa,where the parameter function l1(a) represents thefirst Lyapunov coefficient. If there existcritical values a0 of the parameters, whereC1-1βa0=μa0ω0=0,l1a0=Rec1a0ω0≠0,then, by the linear local (around a0) transformation w=u/l1(a)1/2, (43) yields(44)dudτ=β+iu+suu2+Ou4,with s=±1=sign⁡l1(a0)=sign⁡Re[c1(a0)]. According to [14] (Lemma 3.2, valid for both the supercritical and subcritical case) the terms which are higher than the third order in (44) do not affect the bifurcation behavior of the system near x0, and thus we consider thelocally topologically equivalent system (44), where the O(u4) terms have been dropped. Moreover, if thetransversality condition holds at the critical values of the parameters, that is,C1-2∂β∂aia=a0=1ω0∂μ∂aia=a0≠0,i=1,2,3,then one limit cycle is bifurcated from the equilibrium x0 at a0. By using polar coordinates u=ρeiφ, the equivalent system (44) (without the O(u4) terms) gives dρ/dτ=ρ(β∓ρ2), s=∓1, corresponding to the well-known supercritical (s=-1, β>0) and subcritical (s=1, β<0) cases of the Hopf bifurcation. We call x0(a) a Hopf point ofcodimension 1 (H1 point).In the case wherel1(a0)=0, we proceed to the fifth-order normal form.


#### 3.1.2. Normal Form of the Fifth Order(nr=2)

Here a fifth-ordernormal form is obtained; namely,(45)w˙=λaw+c1aww2+c2aww4+Ow6,with c2(a)w3w-2, the fifth-orderresonant term, due to the term (λ+λ-)h32 included in K32 (we set h32=0). Then by means of the linear time scaling (41) as well as the nonlinear time reparametrization(46)dθ=1-Imc1aωaw2-Imc2aωa-Imc1a2ω2aw4dτ,(45) takes the form of thelocally orbitally equivalent system(47)dwdτ=β+iw+l1aww2+l2aww4+Ow6,l2a=Rec2aωa-βaImc2aωa-Imc1aωal1awith β(a) and l1(a) as in (43) and l2(a), thesecond Lyapunov coefficient. Here if there existcritical values a0 of the parameters, whereC2-1βa0=0,l1a0=0,l2a0=Rec2a0ω0≠0,then the linear local transformation w=u/l2(a)1/4 transforms (47) into(48)dudτ=β+iu+β1uu2+suu4+Ou6,β1a=l1al2a1/2,where s=±1=sign⁡l2(a0)=sign⁡Rec2(a0). For the degenerate case under consideration (Bautin bifurcation), by dropping the Ou6 terms in (48), alocal topological equivalence can be established near x0 ([14], Lemma 8.5). Furthermore, thetransversality condition is expressed by the local smooth invertibility of the map a→(μ(a),l1(a))T at a=a0; namely,C2-2det⁡∂β,l1∂ai,aja=a0=1ω0det⁡∂μ,l1∂ai,aja=a0≠0,i=1,2,j=i+1,…,3,where(49)∂l1∂aia=a0=1ω0Re∂c1∂aia=a0-Imc1a0ω02∂μ∂aia=a0,i=1,2,3.The resulting bifurcation behavior in the parameter space (β,β1) is given below (Section 3.2.1). The equilibrium x0(a) represents a Hopf point ofcodimension 2 (H2 point).If the second Lyapunov coefficient also vanishes ata0 (l2(a0)=0), we proceed to the seventh-order normal form.


#### 3.1.3. Normal Form of the Seventh Order(nr=3)

In this case a seventh-orderresonant term further appears, due to the term (λ+λ-)h43 included in K43 (we set h43=0) and thus the correspondingnormal form becomes(50)w˙=λaw+c1aww2+c2aww4+c3aww6+Ow8,where application of the linear time scaling (41) and the nonlinear time reparametrization(51)dθ=1-Imc1aωaw2-Imc2aωa-Imc1a2ω2aw4-Imc3aωa-2Imc1aImc2aω2a+Imc1a3ω3aw6dτresults in the locally orbitally equivalent system(52)dwdτ=β+iw+l1aww2+l2aww4+l3aww6+Ow8,l3a=Rec3aωa-βaImc3aωa-Imc1aωal2a-Imc2aωal1awith β(a), l1(a), and l2(a) as in (43) and (47), respectively, and l3(a), thethird Lyapunov coefficient. Moreover, the correspondingcritical condition,C3-1βa0=0,l1a0=0,l2a0=0,l3a0=Rec3a0ω0≠0,allows the linear local transformation w=u/l3(a)1/6 to transform (51) as(53)dudτ=β+iu+β13uu2+β2uu4+suu6+Ou8,β13a=l1al3a1/3,β2a=l2al3a1/3,with s=±1=sign⁡l3(a0)=sign⁡Rec3(a0). Alocally topologically equivalent system is also obtained here by dropping the Ou8 terms in (53). Additionally thetransversality condition referring to this degenerate case is expressed by the relationsC3-2adet∂β,l1∂ai,aja=a0=1ω0det∂μ,l1∂ai,aja=a0≠0,i=1,2,j=i+1,…,3,C3-2bdet∂β,l2∂ai,aja=a0=1ω0det∂μ,l2∂ai,aja=a0≠0,i=1,2,j=i+1,…,3,C3-2cdet∂l1,l2∂ai,aja=a0≠0,i=1,2,j=i+1,…,3,where the derivatives of l1 with respect to the parameters are given by (49) and(54)∂l2∂aia=a0=1ω0Re∂c2∂aia=a0-Imc2a0ω02∂μ∂aia=a0-Imc1a0ω0∂l1∂aia=a0,i=1,2,3.The corresponding bifurcation portrait in the parameter space (β,β13,β2) is presented in Section 3.2.2. The equilibrium x0(a) is called a Hopf point ofcodimension 3 (H3 point).Through the same procedure, the vanishing of the derivedLyapunov coefficients in each case leads to higher order normal forms of the system.


### 3.2. Bifurcation Portraits of Degenerate Hopf Bifurcations

#### 3.2.1. Bifurcation of H2 Point

By using polar coordinatesu=ρeiφ, system (48) (without the Ou6 terms) results in(55)dρdτ=∓ρAρ,Aρ=ρ4∓β1ρ2∓β,s=∓1.By denoting Δ=β12±4β for s=∓1, respectively, and investigating the sign of A as a quadratic expression with respect to ρ2, we arrive at zero, one, or two circular limit cycles. More precisely, we obtain the following regions in the (β,β1) plane.Case 1 (s=-1 (Δ=β12+4β)).
(i) RegionR0: Δ<0(β<0) or Δ=0 or β<0, β1<0, Δ>0.
TheBautin pointx0(a) is a (nonlinearly) stable equilibrium (no cycles).
(ii) RegionR1: β>0.
One stable limit cycle with radiusρ2=((β1+Δ)/2)1/2 is bifurcated locally around the unstable equilibrium.
(iii) RegionR2: β<0, β1>0, Δ>0.
We have two homocentric limit cycles with radiusρ1,2=((β1∓Δ)/2)1/2. The outer cycle (ρ=ρ2) is stable, while the inner cycle (ρ=ρ1) is unstable, surrounding the stable equilibrium.For the opposite sign ofs a symmetric bifurcation diagram is generated.Case 2 (s=1 (Δ=β12-4β)).
(i) RegionR0: Δ<0(β>0) or Δ=0 or β>0, β1>0, Δ>0.
TheBautin pointx0(a) is an (nonlinearly) unstable equilibrium (no cycles).
(ii) RegionR1: β<0.
One unstable limit cycle with radiusρ2=((-β1+Δ)/2)1/2 is bifurcated locally around the stable equilibrium.
(iii) RegionR2: β>0, β1<0, Δ>0.
We have two homocentric limit cycles with radiusρ1,2=((-β1∓Δ)/2)1/2. The outer cycle (ρ=ρ2) is unstable, while the inner cycle (ρ=ρ1) is stable, surrounding the unstable equilibrium.


#### 3.2.2. Bifurcation of H3 Point

In this case, by means of the topologically equivalent system (53), we also obtain (55) with A(ρ)=ρ6∓β2ρ4∓β13ρ2∓β, s=∓1. This cubic expression, with respect to ρ2, is reduced to the form A(ρ)=A1(r)=r3+P1r+P2 by means of the substitution ρ2=r±β2/3, where(56)P1=-β223∓β13,P2=∓2β2327-β13β23∓β,s=∓1.As it is well known, the sign of the discriminant Δ1=P13/27+P22/4 determines the kind of the roots of A1. Furthermore, investigation of the sign of A yields the bifurcation regions in the (β,β13,β2) space, which correspond to this degenerate case. Thus by denoting the simple real roots of the cubic A(ρ) as(57)Δ1>0:ρ+2=-P22+Δ13+-P22-Δ13±β23,Δ1=0:ρ02=2-P223±β23,Δ1<0:ρ-12=2-P13cos⁡α3±β23,ρ-2,32=-2-P13cos⁡α3±π3±β23,α=arccos⁡-P22-p3/27,ρ-12>ρ-22>ρ-32,we conclude with zero, one, two, or three circular limit cycles, as follows:(i) RegionR0: Δ1>0, ρ+2<0 or Δ1=0, ρ02<0 or Δ1<0, 0>ρ-12.The pointx0(a) is a nonlinearly stable (unstable) equilibrium for s=-1 (s=1) (no cycles).(ii) RegionR1: Δ1>0, ρ+2>0 or Δ1=0, ρ02>0 or Δ1<0, ρ-12>0>ρ-22>ρ-32.One stable (unstable) limit cycle with radius,(58)ρ1=ρ+2,Δ1>0,ρ02,Δ1=0,ρ-12,Δ1<0,is bifurcated locally around the unstable (stable) equilibrium for s=-1 (s=1).(iii) RegionR2: Δ1<0, ρ-12>ρ-22>0>ρ-32.Here two homocentric limit cycles with radiusρi=ρ-i2, i=1,2, are bifurcated around the stable (unstable) equilibrium for s=-1 (s=1). The outer one (ρ=ρ1) is stable (unstable), while the inner one (ρ=ρ2) is unstable (stable) for s=-1 (s=1).(iv) RegionR3: Δ1<0, ρ-12>ρ-22>ρ-32>0.In this case we have three homocentric limit cycles with radiusρi=ρ-i2, i=1,2,3, surrounding the unstable (stable) equilibrium for s=-1 (s=1). The outer cycle (ρ=ρ1) is stable (unstable) and the middle cycle (ρ=ρ2) is unstable (stable), while the inner one (ρ=ρ3) is stable (unstable) for s=-1 (s=1).


## 3.1. Planar Case(n=2): Poincaré Normal Forms, Codimension nr

In this case, the termsgijm+ with i+j=0,1,2,…, m+=1,2,… in the first equation of (15) and in (16) vanish. Moreover, (17) does not exist. Now, in order to obtain the desired normal form, we first introduce the transformation(35)z=w+h2w,w-+⋯+h2nr+1w,w-,nr∈N,hmw,w-=∑k+l=m1k!l!hklwkw-l,m=2,…,2nr+1,the inversion of which is given by (see [14])(36)w=z-h2z,z--⋯-h2nr+1z,z-+Oz2nr+2.Then by differentiating (36) with respect to t and substituting the right-hand side of (11) for z˙ (with g=g2+g3+⋯), as well as the corresponding conjugate equation for z-˙, by keeping terms up to 2nr+1 order with respect to z, we obtain(37)w˙=λz+g2+⋯+g2nr+1-∑m=12nr-1λz+g2+⋯+g2nr+1-mHm-λ-z-+g-2+⋯+g-2nr+1-mH~m-λzH2nr-λ-z-H~2nr+Oz2nr+2,Hmz,z-=∂hm+1z,z-∂z,H~mz,z-=∂hm+1z,z-∂z-,m=1,…,2nr.Finally substitution of transformation (35) for z (as well as the respective conjugate relation for z-) in (37) yields(38)w˙=λaw+K2+⋯+K2nr+1+Ow2nr+2,Kmw,w-;a=∑k+l=mKklλ,λ-,hkl,gkl,hij,h-ij,gij,g-ijwkw-l,m=2,…,2nr+1,i+j=2,…,m-1,where terms up to 2nr+1 order with respect to w have been evaluated. Now, by setting Kkl, k+l=2,…,2nr+1, equal to zero we define the coefficients hkl of the transformation (35), while the remainingresonant terms of the odd order, ci=Ki+1,i, i=1,…,nr, constitute the desirable normal form of (11). Furthermore, by substituting the obtained expressions of hkl,h-kl involved in ci (k+l=2,…,2i), then by means of computer calculations we arrive at the analytic formulae of the normal coefficients:(39)ci=ciλ,λ-,gi+1,i,gst,g-st,s+t=2,…,2i,i=1,…,nr,which are lengthy for i≥2.


### 3.1.1. Normal Form of the Third Order(nr=1)

In this case, the above procedure results in the third-ordernormal form:(40)w˙=λaw+c1aww2+Ow4,where c1(a)w2w- represents the cubicresonant term, due to the term (λ+λ-)h21 included in K21 (we set h21=0). Then by means of the linear time scaling(41)θ=ωatand the nonlinear time reparametrization(42)dθ=1-Imc1aωaw2dτ, (40) becomes(43)dwdτ=β+iw+l1aww2+Ow4,βa=μaωa,l1a=Rec1aωa-βaImc1aωa,where the parameter function l1(a) represents thefirst Lyapunov coefficient. If there existcritical values a0 of the parameters, whereC1-1βa0=μa0ω0=0,l1a0=Rec1a0ω0≠0,then, by the linear local (around a0) transformation w=u/l1(a)1/2, (43) yields(44)dudτ=β+iu+suu2+Ou4,with s=±1=sign⁡l1(a0)=sign⁡Re[c1(a0)]. According to [14] (Lemma 3.2, valid for both the supercritical and subcritical case) the terms which are higher than the third order in (44) do not affect the bifurcation behavior of the system near x0, and thus we consider thelocally topologically equivalent system (44), where the O(u4) terms have been dropped. Moreover, if thetransversality condition holds at the critical values of the parameters, that is,C1-2∂β∂aia=a0=1ω0∂μ∂aia=a0≠0,i=1,2,3,then one limit cycle is bifurcated from the equilibrium x0 at a0. By using polar coordinates u=ρeiφ, the equivalent system (44) (without the O(u4) terms) gives dρ/dτ=ρ(β∓ρ2), s=∓1, corresponding to the well-known supercritical (s=-1, β>0) and subcritical (s=1, β<0) cases of the Hopf bifurcation. We call x0(a) a Hopf point ofcodimension 1 (H1 point).In the case wherel1(a0)=0, we proceed to the fifth-order normal form.


### 3.1.2. Normal Form of the Fifth Order(nr=2)

Here a fifth-ordernormal form is obtained; namely,(45)w˙=λaw+c1aww2+c2aww4+Ow6,with c2(a)w3w-2, the fifth-orderresonant term, due to the term (λ+λ-)h32 included in K32 (we set h32=0). Then by means of the linear time scaling (41) as well as the nonlinear time reparametrization(46)dθ=1-Imc1aωaw2-Imc2aωa-Imc1a2ω2aw4dτ,(45) takes the form of thelocally orbitally equivalent system(47)dwdτ=β+iw+l1aww2+l2aww4+Ow6,l2a=Rec2aωa-βaImc2aωa-Imc1aωal1awith β(a) and l1(a) as in (43) and l2(a), thesecond Lyapunov coefficient. Here if there existcritical values a0 of the parameters, whereC2-1βa0=0,l1a0=0,l2a0=Rec2a0ω0≠0,then the linear local transformation w=u/l2(a)1/4 transforms (47) into(48)dudτ=β+iu+β1uu2+suu4+Ou6,β1a=l1al2a1/2,where s=±1=sign⁡l2(a0)=sign⁡Rec2(a0). For the degenerate case under consideration (Bautin bifurcation), by dropping the Ou6 terms in (48), alocal topological equivalence can be established near x0 ([14], Lemma 8.5). Furthermore, thetransversality condition is expressed by the local smooth invertibility of the map a→(μ(a),l1(a))T at a=a0; namely,C2-2det⁡∂β,l1∂ai,aja=a0=1ω0det⁡∂μ,l1∂ai,aja=a0≠0,i=1,2,j=i+1,…,3,where(49)∂l1∂aia=a0=1ω0Re∂c1∂aia=a0-Imc1a0ω02∂μ∂aia=a0,i=1,2,3.The resulting bifurcation behavior in the parameter space (β,β1) is given below (Section 3.2.1). The equilibrium x0(a) represents a Hopf point ofcodimension 2 (H2 point).If the second Lyapunov coefficient also vanishes ata0 (l2(a0)=0), we proceed to the seventh-order normal form.


### 3.1.3. Normal Form of the Seventh Order(nr=3)

In this case a seventh-orderresonant term further appears, due to the term (λ+λ-)h43 included in K43 (we set h43=0) and thus the correspondingnormal form becomes(50)w˙=λaw+c1aww2+c2aww4+c3aww6+Ow8,where application of the linear time scaling (41) and the nonlinear time reparametrization(51)dθ=1-Imc1aωaw2-Imc2aωa-Imc1a2ω2aw4-Imc3aωa-2Imc1aImc2aω2a+Imc1a3ω3aw6dτresults in the locally orbitally equivalent system(52)dwdτ=β+iw+l1aww2+l2aww4+l3aww6+Ow8,l3a=Rec3aωa-βaImc3aωa-Imc1aωal2a-Imc2aωal1awith β(a), l1(a), and l2(a) as in (43) and (47), respectively, and l3(a), thethird Lyapunov coefficient. Moreover, the correspondingcritical condition,C3-1βa0=0,l1a0=0,l2a0=0,l3a0=Rec3a0ω0≠0,allows the linear local transformation w=u/l3(a)1/6 to transform (51) as(53)dudτ=β+iu+β13uu2+β2uu4+suu6+Ou8,β13a=l1al3a1/3,β2a=l2al3a1/3,with s=±1=sign⁡l3(a0)=sign⁡Rec3(a0). Alocally topologically equivalent system is also obtained here by dropping the Ou8 terms in (53). Additionally thetransversality condition referring to this degenerate case is expressed by the relationsC3-2adet∂β,l1∂ai,aja=a0=1ω0det∂μ,l1∂ai,aja=a0≠0,i=1,2,j=i+1,…,3,C3-2bdet∂β,l2∂ai,aja=a0=1ω0det∂μ,l2∂ai,aja=a0≠0,i=1,2,j=i+1,…,3,C3-2cdet∂l1,l2∂ai,aja=a0≠0,i=1,2,j=i+1,…,3,where the derivatives of l1 with respect to the parameters are given by (49) and(54)∂l2∂aia=a0=1ω0Re∂c2∂aia=a0-Imc2a0ω02∂μ∂aia=a0-Imc1a0ω0∂l1∂aia=a0,i=1,2,3.The corresponding bifurcation portrait in the parameter space (β,β13,β2) is presented in Section 3.2.2. The equilibrium x0(a) is called a Hopf point ofcodimension 3 (H3 point).Through the same procedure, the vanishing of the derivedLyapunov coefficients in each case leads to higher order normal forms of the system.


## 3.1.1. Normal Form of the Third Order(nr=1)

In this case, the above procedure results in the third-ordernormal form:(40)w˙=λaw+c1aww2+Ow4,where c1(a)w2w- represents the cubicresonant term, due to the term (λ+λ-)h21 included in K21 (we set h21=0). Then by means of the linear time scaling(41)θ=ωatand the nonlinear time reparametrization(42)dθ=1-Imc1aωaw2dτ, (40) becomes(43)dwdτ=β+iw+l1aww2+Ow4,βa=μaωa,l1a=Rec1aωa-βaImc1aωa,where the parameter function l1(a) represents thefirst Lyapunov coefficient. If there existcritical values a0 of the parameters, whereC1-1βa0=μa0ω0=0,l1a0=Rec1a0ω0≠0,then, by the linear local (around a0) transformation w=u/l1(a)1/2, (43) yields(44)dudτ=β+iu+suu2+Ou4,with s=±1=sign⁡l1(a0)=sign⁡Re[c1(a0)]. According to [14] (Lemma 3.2, valid for both the supercritical and subcritical case) the terms which are higher than the third order in (44) do not affect the bifurcation behavior of the system near x0, and thus we consider thelocally topologically equivalent system (44), where the O(u4) terms have been dropped. Moreover, if thetransversality condition holds at the critical values of the parameters, that is,C1-2∂β∂aia=a0=1ω0∂μ∂aia=a0≠0,i=1,2,3,then one limit cycle is bifurcated from the equilibrium x0 at a0. By using polar coordinates u=ρeiφ, the equivalent system (44) (without the O(u4) terms) gives dρ/dτ=ρ(β∓ρ2), s=∓1, corresponding to the well-known supercritical (s=-1, β>0) and subcritical (s=1, β<0) cases of the Hopf bifurcation. We call x0(a) a Hopf point ofcodimension 1 (H1 point).In the case wherel1(a0)=0, we proceed to the fifth-order normal form.


## 3.1.2. Normal Form of the Fifth Order(nr=2)

Here a fifth-ordernormal form is obtained; namely,(45)w˙=λaw+c1aww2+c2aww4+Ow6,with c2(a)w3w-2, the fifth-orderresonant term, due to the term (λ+λ-)h32 included in K32 (we set h32=0). Then by means of the linear time scaling (41) as well as the nonlinear time reparametrization(46)dθ=1-Imc1aωaw2-Imc2aωa-Imc1a2ω2aw4dτ,(45) takes the form of thelocally orbitally equivalent system(47)dwdτ=β+iw+l1aww2+l2aww4+Ow6,l2a=Rec2aωa-βaImc2aωa-Imc1aωal1awith β(a) and l1(a) as in (43) and l2(a), thesecond Lyapunov coefficient. Here if there existcritical values a0 of the parameters, whereC2-1βa0=0,l1a0=0,l2a0=Rec2a0ω0≠0,then the linear local transformation w=u/l2(a)1/4 transforms (47) into(48)dudτ=β+iu+β1uu2+suu4+Ou6,β1a=l1al2a1/2,where s=±1=sign⁡l2(a0)=sign⁡Rec2(a0). For the degenerate case under consideration (Bautin bifurcation), by dropping the Ou6 terms in (48), alocal topological equivalence can be established near x0 ([14], Lemma 8.5). Furthermore, thetransversality condition is expressed by the local smooth invertibility of the map a→(μ(a),l1(a))T at a=a0; namely,C2-2det⁡∂β,l1∂ai,aja=a0=1ω0det⁡∂μ,l1∂ai,aja=a0≠0,i=1,2,j=i+1,…,3,where(49)∂l1∂aia=a0=1ω0Re∂c1∂aia=a0-Imc1a0ω02∂μ∂aia=a0,i=1,2,3.The resulting bifurcation behavior in the parameter space (β,β1) is given below (Section 3.2.1). The equilibrium x0(a) represents a Hopf point ofcodimension 2 (H2 point).If the second Lyapunov coefficient also vanishes ata0 (l2(a0)=0), we proceed to the seventh-order normal form.


## 3.1.3. Normal Form of the Seventh Order(nr=3)

In this case a seventh-orderresonant term further appears, due to the term (λ+λ-)h43 included in K43 (we set h43=0) and thus the correspondingnormal form becomes(50)w˙=λaw+c1aww2+c2aww4+c3aww6+Ow8,where application of the linear time scaling (41) and the nonlinear time reparametrization(51)dθ=1-Imc1aωaw2-Imc2aωa-Imc1a2ω2aw4-Imc3aωa-2Imc1aImc2aω2a+Imc1a3ω3aw6dτresults in the locally orbitally equivalent system(52)dwdτ=β+iw+l1aww2+l2aww4+l3aww6+Ow8,l3a=Rec3aωa-βaImc3aωa-Imc1aωal2a-Imc2aωal1awith β(a), l1(a), and l2(a) as in (43) and (47), respectively, and l3(a), thethird Lyapunov coefficient. Moreover, the correspondingcritical condition,C3-1βa0=0,l1a0=0,l2a0=0,l3a0=Rec3a0ω0≠0,allows the linear local transformation w=u/l3(a)1/6 to transform (51) as(53)dudτ=β+iu+β13uu2+β2uu4+suu6+Ou8,β13a=l1al3a1/3,β2a=l2al3a1/3,with s=±1=sign⁡l3(a0)=sign⁡Rec3(a0). Alocally topologically equivalent system is also obtained here by dropping the Ou8 terms in (53). Additionally thetransversality condition referring to this degenerate case is expressed by the relationsC3-2adet∂β,l1∂ai,aja=a0=1ω0det∂μ,l1∂ai,aja=a0≠0,i=1,2,j=i+1,…,3,C3-2bdet∂β,l2∂ai,aja=a0=1ω0det∂μ,l2∂ai,aja=a0≠0,i=1,2,j=i+1,…,3,C3-2cdet∂l1,l2∂ai,aja=a0≠0,i=1,2,j=i+1,…,3,where the derivatives of l1 with respect to the parameters are given by (49) and(54)∂l2∂aia=a0=1ω0Re∂c2∂aia=a0-Imc2a0ω02∂μ∂aia=a0-Imc1a0ω0∂l1∂aia=a0,i=1,2,3.The corresponding bifurcation portrait in the parameter space (β,β13,β2) is presented in Section 3.2.2. The equilibrium x0(a) is called a Hopf point ofcodimension 3 (H3 point).Through the same procedure, the vanishing of the derivedLyapunov coefficients in each case leads to higher order normal forms of the system.


## 3.2. Bifurcation Portraits of Degenerate Hopf Bifurcations

### 3.2.1. Bifurcation of H2 Point

By using polar coordinatesu=ρeiφ, system (48) (without the Ou6 terms) results in(55)dρdτ=∓ρAρ,Aρ=ρ4∓β1ρ2∓β,s=∓1.By denoting Δ=β12±4β for s=∓1, respectively, and investigating the sign of A as a quadratic expression with respect to ρ2, we arrive at zero, one, or two circular limit cycles. More precisely, we obtain the following regions in the (β,β1) plane.Case 1 (s=-1 (Δ=β12+4β)).
(i) RegionR0: Δ<0(β<0) or Δ=0 or β<0, β1<0, Δ>0.
TheBautin pointx0(a) is a (nonlinearly) stable equilibrium (no cycles).
(ii) RegionR1: β>0.
One stable limit cycle with radiusρ2=((β1+Δ)/2)1/2 is bifurcated locally around the unstable equilibrium.
(iii) RegionR2: β<0, β1>0, Δ>0.
We have two homocentric limit cycles with radiusρ1,2=((β1∓Δ)/2)1/2. The outer cycle (ρ=ρ2) is stable, while the inner cycle (ρ=ρ1) is unstable, surrounding the stable equilibrium.For the opposite sign ofs a symmetric bifurcation diagram is generated.Case 2 (s=1 (Δ=β12-4β)).
(i) RegionR0: Δ<0(β>0) or Δ=0 or β>0, β1>0, Δ>0.
TheBautin pointx0(a) is an (nonlinearly) unstable equilibrium (no cycles).
(ii) RegionR1: β<0.
One unstable limit cycle with radiusρ2=((-β1+Δ)/2)1/2 is bifurcated locally around the stable equilibrium.
(iii) RegionR2: β>0, β1<0, Δ>0.
We have two homocentric limit cycles with radiusρ1,2=((-β1∓Δ)/2)1/2. The outer cycle (ρ=ρ2) is unstable, while the inner cycle (ρ=ρ1) is stable, surrounding the unstable equilibrium.


### 3.2.2. Bifurcation of H3 Point

In this case, by means of the topologically equivalent system (53), we also obtain (55) with A(ρ)=ρ6∓β2ρ4∓β13ρ2∓β, s=∓1. This cubic expression, with respect to ρ2, is reduced to the form A(ρ)=A1(r)=r3+P1r+P2 by means of the substitution ρ2=r±β2/3, where(56)P1=-β223∓β13,P2=∓2β2327-β13β23∓β,s=∓1.As it is well known, the sign of the discriminant Δ1=P13/27+P22/4 determines the kind of the roots of A1. Furthermore, investigation of the sign of A yields the bifurcation regions in the (β,β13,β2) space, which correspond to this degenerate case. Thus by denoting the simple real roots of the cubic A(ρ) as(57)Δ1>0:ρ+2=-P22+Δ13+-P22-Δ13±β23,Δ1=0:ρ02=2-P223±β23,Δ1<0:ρ-12=2-P13cos⁡α3±β23,ρ-2,32=-2-P13cos⁡α3±π3±β23,α=arccos⁡-P22-p3/27,ρ-12>ρ-22>ρ-32,we conclude with zero, one, two, or three circular limit cycles, as follows:(i) RegionR0: Δ1>0, ρ+2<0 or Δ1=0, ρ02<0 or Δ1<0, 0>ρ-12.The pointx0(a) is a nonlinearly stable (unstable) equilibrium for s=-1 (s=1) (no cycles).(ii) RegionR1: Δ1>0, ρ+2>0 or Δ1=0, ρ02>0 or Δ1<0, ρ-12>0>ρ-22>ρ-32.One stable (unstable) limit cycle with radius,(58)ρ1=ρ+2,Δ1>0,ρ02,Δ1=0,ρ-12,Δ1<0,is bifurcated locally around the unstable (stable) equilibrium for s=-1 (s=1).(iii) RegionR2: Δ1<0, ρ-12>ρ-22>0>ρ-32.Here two homocentric limit cycles with radiusρi=ρ-i2, i=1,2, are bifurcated around the stable (unstable) equilibrium for s=-1 (s=1). The outer one (ρ=ρ1) is stable (unstable), while the inner one (ρ=ρ2) is unstable (stable) for s=-1 (s=1).(iv) RegionR3: Δ1<0, ρ-12>ρ-22>ρ-32>0.In this case we have three homocentric limit cycles with radiusρi=ρ-i2, i=1,2,3, surrounding the unstable (stable) equilibrium for s=-1 (s=1). The outer cycle (ρ=ρ1) is stable (unstable) and the middle cycle (ρ=ρ2) is unstable (stable), while the inner one (ρ=ρ3) is stable (unstable) for s=-1 (s=1).


## 3.2.1. Bifurcation of H2 Point

By using polar coordinatesu=ρeiφ, system (48) (without the Ou6 terms) results in(55)dρdτ=∓ρAρ,Aρ=ρ4∓β1ρ2∓β,s=∓1.By denoting Δ=β12±4β for s=∓1, respectively, and investigating the sign of A as a quadratic expression with respect to ρ2, we arrive at zero, one, or two circular limit cycles. More precisely, we obtain the following regions in the (β,β1) plane.Case 1 (s=-1 (Δ=β12+4β)).
(i) RegionR0: Δ<0(β<0) or Δ=0 or β<0, β1<0, Δ>0.
TheBautin pointx0(a) is a (nonlinearly) stable equilibrium (no cycles).
(ii) RegionR1: β>0.
One stable limit cycle with radiusρ2=((β1+Δ)/2)1/2 is bifurcated locally around the unstable equilibrium.
(iii) RegionR2: β<0, β1>0, Δ>0.
We have two homocentric limit cycles with radiusρ1,2=((β1∓Δ)/2)1/2. The outer cycle (ρ=ρ2) is stable, while the inner cycle (ρ=ρ1) is unstable, surrounding the stable equilibrium.For the opposite sign ofs a symmetric bifurcation diagram is generated.Case 2 (s=1 (Δ=β12-4β)).
(i) RegionR0: Δ<0(β>0) or Δ=0 or β>0, β1>0, Δ>0.
TheBautin pointx0(a) is an (nonlinearly) unstable equilibrium (no cycles).
(ii) RegionR1: β<0.
One unstable limit cycle with radiusρ2=((-β1+Δ)/2)1/2 is bifurcated locally around the stable equilibrium.
(iii) RegionR2: β>0, β1<0, Δ>0.
We have two homocentric limit cycles with radiusρ1,2=((-β1∓Δ)/2)1/2. The outer cycle (ρ=ρ2) is unstable, while the inner cycle (ρ=ρ1) is stable, surrounding the unstable equilibrium.


## 3.2.2. Bifurcation of H3 Point

In this case, by means of the topologically equivalent system (53), we also obtain (55) with A(ρ)=ρ6∓β2ρ4∓β13ρ2∓β, s=∓1. This cubic expression, with respect to ρ2, is reduced to the form A(ρ)=A1(r)=r3+P1r+P2 by means of the substitution ρ2=r±β2/3, where(56)P1=-β223∓β13,P2=∓2β2327-β13β23∓β,s=∓1.As it is well known, the sign of the discriminant Δ1=P13/27+P22/4 determines the kind of the roots of A1. Furthermore, investigation of the sign of A yields the bifurcation regions in the (β,β13,β2) space, which correspond to this degenerate case. Thus by denoting the simple real roots of the cubic A(ρ) as(57)Δ1>0:ρ+2=-P22+Δ13+-P22-Δ13±β23,Δ1=0:ρ02=2-P223±β23,Δ1<0:ρ-12=2-P13cos⁡α3±β23,ρ-2,32=-2-P13cos⁡α3±π3±β23,α=arccos⁡-P22-p3/27,ρ-12>ρ-22>ρ-32,we conclude with zero, one, two, or three circular limit cycles, as follows:(i) RegionR0: Δ1>0, ρ+2<0 or Δ1=0, ρ02<0 or Δ1<0, 0>ρ-12.The pointx0(a) is a nonlinearly stable (unstable) equilibrium for s=-1 (s=1) (no cycles).(ii) RegionR1: Δ1>0, ρ+2>0 or Δ1=0, ρ02>0 or Δ1<0, ρ-12>0>ρ-22>ρ-32.One stable (unstable) limit cycle with radius,(58)ρ1=ρ+2,Δ1>0,ρ02,Δ1=0,ρ-12,Δ1<0,is bifurcated locally around the unstable (stable) equilibrium for s=-1 (s=1).(iii) RegionR2: Δ1<0, ρ-12>ρ-22>0>ρ-32.Here two homocentric limit cycles with radiusρi=ρ-i2, i=1,2, are bifurcated around the stable (unstable) equilibrium for s=-1 (s=1). The outer one (ρ=ρ1) is stable (unstable), while the inner one (ρ=ρ2) is unstable (stable) for s=-1 (s=1).(iv) RegionR3: Δ1<0, ρ-12>ρ-22>ρ-32>0.In this case we have three homocentric limit cycles with radiusρi=ρ-i2, i=1,2,3, surrounding the unstable (stable) equilibrium for s=-1 (s=1). The outer cycle (ρ=ρ1) is stable (unstable) and the middle cycle (ρ=ρ2) is unstable (stable), while the inner one (ρ=ρ3) is stable (unstable) for s=-1 (s=1).


## 4. Application in a ModifiedLorenz System

As an example we apply the above analysis and the respective Mathematica software in the system(59)x˙=-σx+σy,y˙=σ-cx+cy+b2xz,z˙=-bz+xy,with (σ,c,b2)∈R3, which can be considered as a modified more general form of theLorenz system [9]. By investigating the stability of the equilibrium (60)E=-σbb2,-σbb2,-σb2,firstly we obtain the characteristic equation(61)λ3+β2λ2+β1λ+β0=0,β2=σ-c+b,β1=2σ-cb,β0=2σ2b.Then by combining the Routh-Hurwitz stability conditions, (62)β0>0,β1>0,β1β2-β0>0,with the relations resulting from (61) for λ=iω0, we arrive at the following necessary and sufficient Hopf conditions related to the equilibrium E:(63)c>0,σc>12orc<0,13<σc<12,(64)b=bcr=c3σ-c2σ-c,ω=ω0=c3σ-c,while the (critical) real root of the characteristic is(65)λ3=-2σ22σ-c<0.We should additionally note that, in the (σ,c) regions defined by (63), E is stable for b>bcr. Thus in the three-dimensional parameter space a1,a2,a3=σ,c,b, the essential feature of the Hopf bifurcation, that is, the appearance of purely imaginary eigenvalues together with a negative real one, occurs at the critical surface a3=a3cr(a1,a2).By following now the steps presented in Section2.3, by means of the corresponding constructed symbolic algorithms, firstly we arrive at the portrait of the codimension 1 Hopf bifurcation with respect to E equilibrium path. Both the (σ,c) regions (b=bcr), leading to supercritical (negative first Lyapunov coefficient) and subcritical (positive first Lyapunov coefficient) bifurcations, are presented (Figure 1). Then we specify the points where l1(a0) becomes zero, all corresponding to the codimension 2 case (Figure 2), since for all the obtained respective couples (σ,c) we find that the critical second Lyapunov coefficient l2(a0) is different from zero. We note that for every positive value of σ there exist two degenerate Hopf points in the (σ,c) plane. We should also note that thetransversality conditions C1-2 and C2-2 hold true for allH1 andH2 points, respectively.Figure 1
Sign of the first Lyapunov coefficient atb=bcr (blue for negative and red for positive values) in the parameter plane σ,c for system (59) (b2=-1).Figure 2
H2 points in the parameter plane σ,c for system (59) (b2=-1).Furthermore, the codimension 2 bifurcation portraits (see the corresponding regions in Section3.2.1) are extracted for allH2 points which are presented in Figure 2. In Figures 3 and 4 two portraits are shown for specific couples (σ,c) of codimension 2, corresponding to negative and positive sign of l2(a0), respectively. We should mention here that Region R2 does not appear for anyH2 point as b varies in the considered interval.Figure 3
Bifurcation portrait ofH2 point (s=-1): σ,c,b2=30,4.1861,-1. Region 0 (green), Region 1 (purple) (red dot marks bcr).Figure 4
Bifurcation portrait ofH2 point (s=1): σ,c,b2=2.2,2.4871,-1. Region 0 (green), Region 1 (purple) (red dot marks bcr).Finally we verify the bifurcation results by determining numerically two cycles emerging in theb-equilibrium path of the twoH2 points related to Figures 3 and 4, respectively (Region R1). More specifically, by use of the orthogonal collocation on finite elements method (which does not use numerical integration), we obtain one stable limit cycle surrounding the unstable equilibrium (sign⁡l2(a0)=-1) (Figure 5) and an unstable one bifurcating around the stable equilibrium (sign⁡l2(a0)=1) (Figure 6). The same cycles are generated by the variable-step, variable-order Adams-Bashforth-Moulton predictor-corrector method of orders 1 to 12, which is the standard Matlab routine used to solve nonstiff ODEs, noted as “ode113” in the graphs illustrated in Figures 5 and 6.Figure 5
Codimension 2 (see (48)), s=-1, Region R1. Stable cycle and behavior of inner orbits for σ,c,b,b2=30,4.1861,5.29,-1 (bcr=6.43618) (pink marker marks the beginning). Unstable equilibrium (black marker): 12.5976,12.5976,30. ode113: Adams-Bashforth-Moulton PECE solver of orders 1 to 12 (Matlab).Figure 6
Codimension 2 (see (48)), s=1, Region R1. Unstable cycle and behavior of inner orbits for σ,c,b,b2=2.2,2.4871,5.47,-1 (bcr=5.3476) (pink marker marks the beginning). Stable equilibrium (black marker): 3.469,3.469,2.2. ode113: Adams-Bashforth-Moulton PECE solver of orders 1 to 12 (Matlab).


## 5. Conclusion

The computation of the bifurcation features of a dynamical system is important as regards the behavior of the phase variables under the change of the values of the various parameters involved. The new symbolic representation presented herein and the analytical procedure that followed are proved excellent at handling multiparameter systems and manipulating the resulting lengthy formulae concerning the center manifold coefficients, the coefficients of the restricted equation, the Lyapunov coefficients, and all respective derivatives involved in thetransversality conditions of the bifurcations.By applying the procedure described in Section2 to a modified Lorenz system and using the symbolic computational software Mathematica 7, we evaluated qualitative and quantitative data with respect to the bifurcation of the system very fast. In addition we verified our results by obtaining limit cycles according to theory corresponding to parameter points of codimension 2.Thus under certain improvements and modifications the symbolic algorithms related to the analysis can be applied to any multidimensional and multiparameter system and evaluate a massive amount of data by saving computational time and memory space.



---

*Source: 7658364-2016-06-30.xml*