International Journal of Mathematics and Mathematical Sciences
Volume 2016 (2016), Article ID 7658364, 12 pages
http://dx.doi.org/10.1155/2016/7658364

Research Article

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

Michail P. Markakis and Panagiotis S. Douris

Department of Electrical & Computer Engineering, University of Patras, 26504 Patras, Greece

Received 18 December 2015; Revised 18 May 2016; Accepted 25 May 2016

Academic Editor: Harvinder S. Sidhu

Copyright © 2016 Michail P. Markakis and Panagiotis S. Douris. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The restriction of an -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.

1. Introduction

In the international literature, we encounter a large number of -dimensional nonlinear dynamical systems () 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 to -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 () ([14], Sections 5.4.1, 8.7.1, and 8.7.3) leads to the computation of the critical polynomial coefficients of 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 an -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 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 (Sections 2, 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 of . 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 the resonant 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, the transversality 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 -Dimensional Coordinate Space

Consider a smooth continuous-time three-parameter system with smooth dependence on the parameters:

with , , and with . If is the equilibrium path of (1) and , the Jacobian matrix evaluated at this equilibrium path, then, for small perturbations , expansion of (1) yields

The smooth vector function represents the nonlinear terms of the right-hand side of (2); that is, and

where

with . It is well known that the Hopf bifurcation theorem (see [17], Section 11.2) implies that if, at a critical triplet , has a pair of pure imaginary eigenvalues , , while the real part of the remaining eigenvalues is negative, then an oscillatory type of instability occurs close to , leading to a family of limit cycles. Moreover we conclude that, in a small region around , has a pair of complex conjugate eigenvalues and with

We also note that the classic Hopf theory additionally demands that, at , cross the imaginary axis with nonzero speed (transversality condition).

Let us now consider the complex eigenvectors and of and , respectively, having the properties

and satisfying the normalization condition

By the above normalization it can also be easily proved that

Let denote the two-dimensional parameter-dependent real eigenspace (corresponding to ), spanned by , which becomes critical at , and the ()-dimensional () real eigenspace corresponding to all eigenvalues of other than . Then, by use of the following Lemma ([14], Lemma ) we get the following.

Lemma 1. Consider if and only if .

Consider ; we introduce a complex variable and decompose any as

where and . Taking into account (7) and (8), the last equation yields

Thus taking into account (2), as well as (6) and (9), differentiation of (10) yields the reduced form of system (2) in the -dimensional coordinate space :

where and

Systems (2) and (11)-(12) are in fact dimensionally equivalent, due to the two real constraints imposed on by Lemma 1.

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

The complex function , given by (13), has the following form:

where

with , , , .

Thus by using (14), (11) and (12) become

where

with as in (15), , , , .

2.3. Computation of the Center Manifold

Regarding now the center manifold , we adopt the representation

with , , and . Since we want to investigate the system under consideration up to codimension 3 bifurcation of the equilibrium, by substituting in (16) for the center manifold expression (19) and keeping terms up to seventh order with respect to , we arrive at

where

with . By means of (20) we conclude that we need to evaluate the center manifold coefficients up to the sixth order (, ). Thus by writing the invariance relation for ,

where

we substitute (16) and (17) in (22), where the first equation of (19) is further substituted for in and , , (see (15) and (18)), and keep terms up to the sixth order with respect to . 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, , and for , , and , respectively, we equate the terms of the same order, from the second up to the sixth and obtain five equations, namely, , .(2) In the obtained equations , by substituting the second equation of (19) for , , as well as the associated expressions giving their derivatives with respect to and , for and , , respectively, we conclude with new explicit forms of the above equations, namely, , .(3) In each equation , , by considering the coefficients of the terms, , and taking into account the following relations,

 based on the first equation of (6), with , we solve for :

 Consider . In particular in the critical case , , for , since the coefficient matrix of the corresponding system for becomes with determinant equal to zero, the solution is obtained in combination with the additional condition (see Section 2.4).(4) By substituting the above derived expressions of , , in (20) and keeping terms up to seventh order with respect to , we obtain the coefficients , , of the two-dimensional equation in the plane, representing the restriction of the system on the center manifold; namely,

(5) Finally, by considering the analytic expressions (39) of the resonant terms corresponding to the respective normal forms (, or , depending on the considered bifurcation codimension; see Section 3.1 below), we substitute , , and with the coefficients , , and , 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 corresponding transversality conditions (see , , and (, , and ) in Sections 3.1.1, 3.1.2, and 3.1.3, resp.).

2.4. Derivation of the Critical Center Manifold Coefficients ,

Due to the form of the equations resulting in the vector coefficients , namely,

or

for , in the critical case , , the determinant becomes zero. Here by means of the constraint

at first it can be easily proved that and hence , which, as it can further be proved, is a sufficient and necessary condition for the solvability of the corresponding system; namely,

The expressions of the solutions for the critical coefficients , also satisfy constraint (29). Moreover as regards the expressions of the critical derivatives of , with respect to , involved in the transversality 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

where , are minor determinants of the matrix:

Here, we also provide the expressions of the vectors and :

with , , minor determinants of and

with , , minor determinants of .

3. Normal Forms and Degenerate Bifurcations

3.1. Planar Case : Poincaré Normal Forms, Codimension

In this case, the terms with , 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

the inversion of which is given by (see [14])

Then by differentiating (36) with respect to and substituting the right-hand side of (11) for (with ), as well as the corresponding conjugate equation for , by keeping terms up to order with respect to , we obtain

Finally substitution of transformation (35) for (as well as the respective conjugate relation for ) in (37) yields

where terms up to order with respect to have been evaluated. Now, by setting , , equal to zero we define the coefficients of the transformation (35), while the remaining resonant terms of the odd order, , , constitute the desirable normal form of (11). Furthermore, by substituting the obtained expressions of involved in (), then by means of computer calculations we arrive at the analytic formulae of the normal coefficients:

which are lengthy for .

3.1.1. Normal Form of the Third Order

In this case, the above procedure results in the third-order normal form:

where represents the cubic resonant term, due to the term included in (we set ). Then by means of the linear time scaling

and the nonlinear time reparametrization

(40) becomes

where the parameter function represents the first Lyapunov coefficient. If there exist critical values of the parameters, where

then, by the linear local (around ) transformation , (43) yields

with . 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 , and thus we consider the locally topologically equivalent system (44), where the terms have been dropped. Moreover, if the transversality condition holds at the critical values of the parameters, that is,

then one limit cycle is bifurcated from the equilibrium at . By using polar coordinates , the equivalent system (44) (without the terms) gives , , corresponding to the well-known supercritical (, ) and subcritical (, ) cases of the Hopf bifurcation. We call a Hopf point of codimension 1 (H1 point).

In the case where , we proceed to the fifth-order normal form.

3.1.2. Normal Form of the Fifth Order

Here a fifth-order normal form is obtained; namely,

with , the fifth-order resonant term, due to the term included in (we set ). Then by means of the linear time scaling (41) as well as the nonlinear time reparametrization

(45) takes the form of the locally orbitally equivalent system

with and as in (43) and , the second Lyapunov coefficient. Here if there exist critical values of the parameters, where

then the linear local transformation transforms (47) into

where . For the degenerate case under consideration (Bautin bifurcation), by dropping the terms in (48), a local topological equivalence can be established near ([14], Lemma 8.5). Furthermore, the transversality condition is expressed by the local smooth invertibility of the map at ; namely,

where

The resulting bifurcation behavior in the parameter space is given below (Section 3.2.1). The equilibrium represents a Hopf point of codimension 2 (H2 point).

If the second Lyapunov coefficient also vanishes at (), we proceed to the seventh-order normal form.

3.1.3. Normal Form of the Seventh Order

In this case a seventh-order resonant term further appears, due to the term included in (we set ) and thus the corresponding normal form becomes

where application of the linear time scaling (41) and the nonlinear time reparametrization

results in the locally orbitally equivalent system

with , , and as in (43) and (47), respectively, and , the third Lyapunov coefficient. Moreover, the corresponding critical condition,

allows the linear local transformation to transform (51) as

with . A locally topologically equivalent system is also obtained here by dropping the terms in (53). Additionally the transversality condition referring to this degenerate case is expressed by the relations

where the derivatives of with respect to the parameters are given by (49) and

The corresponding bifurcation portrait in the parameter space is presented in Section 3.2.2. The equilibrium is called a Hopf point of codimension 3 (H3 point).

Through the same procedure, the vanishing of the derived Lyapunov 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 coordinates , system (48) (without the terms) results in

By denoting for , respectively, and investigating the sign of as a quadratic expression with respect to , we arrive at zero, one, or two circular limit cycles. More precisely, we obtain the following regions in the plane.

Case 1 ( ()). (i) Region : or or , , .
The Bautin point is a (nonlinearly) stable equilibrium (no cycles).
(ii) Region : .
One stable limit cycle with radius is bifurcated locally around the unstable equilibrium.
(iii) Region : , , .
We have two homocentric limit cycles with radius . The outer cycle () is stable, while the inner cycle () is unstable, surrounding the stable equilibrium.

For the opposite sign of a symmetric bifurcation diagram is generated.

Case 2 ( ()). (i) Region : or or , , .
The Bautin point is an (nonlinearly) unstable equilibrium (no cycles).
(ii) Region : .
One unstable limit cycle with radius is bifurcated locally around the stable equilibrium.
(iii) Region : , , .
We have two homocentric limit cycles with radius . The outer cycle () is unstable, while the inner cycle () 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 , . This cubic expression, with respect to , is reduced to the form by means of the substitution , where

As it is well known, the sign of the discriminant determines the kind of the roots of . Furthermore, investigation of the sign of yields the bifurcation regions in the space, which correspond to this degenerate case. Thus by denoting the simple real roots of the cubic as

we conclude with zero, one, two, or three circular limit cycles, as follows:

(i) Region : , or , or , .

The point is a nonlinearly stable (unstable) equilibrium for () (no cycles).

(ii) Region : , or , or , .

One stable (unstable) limit cycle with radius,

is bifurcated locally around the unstable (stable) equilibrium for ().

(iii) Region : , .

Here two homocentric limit cycles with radius , , are bifurcated around the stable (unstable) equilibrium for (). The outer one () is stable (unstable), while the inner one () is unstable (stable) for ().

(iv) Region : , .

In this case we have three homocentric limit cycles with radius , , surrounding the unstable (stable) equilibrium for (). The outer cycle () is stable (unstable) and the middle cycle () is unstable (stable), while the inner one () is stable (unstable) for ().

4. Application in a Modified Lorenz System

As an example we apply the above analysis and the respective Mathematica software in the system

with , which can be considered as a modified more general form of the Lorenz system [9]. By investigating the stability of the equilibrium

firstly we obtain the characteristic equation

Then by combining the Routh-Hurwitz stability conditions,

with the relations resulting from (61) for , we arrive at the following necessary and sufficient Hopf conditions related to the equilibrium :

while the (critical) real root of the characteristic is

We should additionally note that, in the regions defined by (63), is stable for . Thus in the three-dimensional parameter space , 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 .

By following now the steps presented in Section 2.3, by means of the corresponding constructed symbolic algorithms, firstly we arrive at the portrait of the codimension 1 Hopf bifurcation with respect to equilibrium path. Both the regions (), leading to supercritical (negative first Lyapunov coefficient) and subcritical (positive first Lyapunov coefficient) bifurcations, are presented (Figure 1). Then we specify the points where becomes zero, all corresponding to the codimension 2 case (Figure 2), since for all the obtained respective couples we find that the critical second Lyapunov coefficient is different from zero. We note that for every positive value of there exist two degenerate Hopf points in the plane. We should also note that the transversality conditions and hold true for all H1 and H2 points, respectively.

Figure 1: Sign of the first Lyapunov coefficient at (blue for negative and red for positive values) in the parameter plane for system (59) ().

Figure 2: H2 points in the parameter plane for system (59) ().

Furthermore, the codimension 2 bifurcation portraits (see the corresponding regions in Section 3.2.1) are extracted for all H2 points which are presented in Figure 2. In Figures 3 and 4 two portraits are shown for specific couples of codimension 2, corresponding to negative and positive sign of , respectively. We should mention here that Region does not appear for any H2 point as varies in the considered interval.

Figure 3: Bifurcation portrait of H2 point (): . Region 0 (green), Region 1 (purple) (red dot marks ).

Figure 4: Bifurcation portrait of H2 point (): . Region 0 (green), Region 1 (purple) (red dot marks ).

Finally we verify the bifurcation results by determining numerically two cycles emerging in the equilibrium path of the two H2 points related to Figures 3 and 4, respectively (Region ). 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 () (Figure 5) and an unstable one bifurcating around the stable equilibrium () (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)), , Region . Stable cycle and behavior of inner orbits for () (pink marker marks the beginning). Unstable equilibrium (black marker): . ode113: Adams-Bashforth-Moulton PECE solver of orders 1 to 12 (Matlab).

Figure 6: Codimension 2 (see (48)), , Region . Unstable cycle and behavior of inner orbits for () (pink marker marks the beginning). Stable equilibrium (black marker): . 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 the transversality conditions of the bifurcations.

By applying the procedure described in Section 2 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.

Competing Interests

The authors declare that they have no competing interests.

References

E. Avila-Vales, N. Chan-Chi, G. E. Garcia-Almeida, and C. Vargas-De-León, “Stability and Hopf bifurcation in a delayed viral infection model with mitosis transmission,” Applied Mathematics and Computation, vol. 259, pp. 293–312, 2015.

W. Wang and S. Ruan, “Bifurcation in an epidemic model with constant removal rate of the infectives,” Journal of Mathematical Analysis and Applications, vol. 291, no. 2, pp. 775–793, 2004.

Z. Feng, Z. Qiu, R. Liu, and D. L. DeAngelis, “Dynamics of a plant-herbivore-predator system with plant-toxicity,” Mathematical Biosciences, vol. 229, no. 2, pp. 190–204, 2011.

Y. Pei, X. Ji, and C. Li, “Pest regulation by means of continuous and impulsive nonlinear controls,” Mathematical and Computer Modelling, vol. 51, no. 5-6, pp. 810–822, 2010.

M. Wang, H. Xu, and Z. Sheng, “A novel three-dimensional energy resources chaotic system with government regulation and empirical research,” Journal of Information and Computational Science, vol. 10, no. 15, pp. 4879–4892, 2013.

M. Sun, L. Tian, and J. Yin, “Hopf bifurcation analysis of the energy resource chaotic system,” International Journal of Nonlinear Science, vol. 1, no. 1, pp. 49–53, 2006.

X. Liu, W. Cai, J. Lu, and Y. Wang, “Stability and Hopf bifurcation for a business cycle model with expectation and delay,” Communications in Nonlinear Science and Numerical Simulation, vol. 25, no. 1–3, pp. 149–161, 2015.

T. Asada, V. Kalantonis, M. Markakis, and P. Markellos, “Analytical expressions of periodic disequilibrium fluctuations generated by Hopf bifurcations in economic dynamics,” Applied Mathematics and Computation, vol. 218, no. 13, pp. 7066–7077, 2012.

E. N. Lorenz, “Deterministic nonperiodic flow,” Journal of the Atmospheric Sciences, vol. 20, pp. 130–141, 1963.

L. F. Mello and S. F. Coelho, “Degenerate Hopf bifurcations in the Lü system,” Physics Letters. A, vol. 373, no. 12-13, pp. 1116–1120, 2009.

Y. Yu and S. Zhang, “Hopf bifurcation analysis of the Lü system,” Chaos, Solitons and Fractals, vol. 21, no. 5, pp. 1215–1220, 2004.

J. Lü, T. Zhou, G. Chen, and S. Zhang, “Local bifurcations of the Chen system,” International Journal of Bifurcation and Chaos, vol. 12, no. 10, pp. 2257–2270, 2002.

K. Huang and Q. Yang, “Stability and Hopf bifurcation analysis of a new system,” Chaos, Solitons and Fractals, vol. 39, no. 2, pp. 567–578, 2009.

Y. A. Kuznetsov, Elements of Applied Bifurcation Theory, Springer, New York, NY, USA, 3rd edition, 2004.

J. Sotomayor, L. F. Mello, and D. de Carvalho Braga, “Lyapunov coefficients for degenerate Hopf bifurcations,” http://arxiv.org/abs/0709.3949

J. Sotomayor, L. F. Mello, and D. D. C. Braga, “Bifurcation analysis of the Watt governor system,” Computational and Applied Mathematics, vol. 26, no. 1, pp. 19–44, 2007.

J. K. Hale and H. Koçak, Dynamics and Bifurcations, Texts in Applied Mathematics, vol. 3, Springer, New York, NY, USA, 1991.