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) yieldsThe smooth vector function represents the nonlinear terms of the right-hand side of (2); that is, andwherewith . 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 withWe 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 propertiesand satisfying the normalization conditionBy the above normalization it can also be easily proved thatLet 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 aswhere and . Taking into account (7) and (8), the last equation yieldsThus 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 andSystems (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:wherewith , , , .
Thus by using (14), (11) and (12) becomewherewith as in (15), , , , .
2.3. Computation of the Center Manifold
Regarding now the center manifold , we adopt the representationwith , , 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 atwherewith . 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 ,wherewe 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,orfor , in the critical case , , the determinant becomes zero. Here by means of the constraintat 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 writewhere , are minor determinants of the matrix:Here, we also provide the expressions of the vectors and :with , , minor determinants of andwith , , 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 transformationthe 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 obtainFinally substitution of transformation (35) for (as well as the respective conjugate relation for ) in (37) yieldswhere 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 scalingand the nonlinear time reparametrization (40) becomeswhere the parameter function represents the first Lyapunov coefficient. If there exist critical values of the parameters, wherethen, by the linear local (around ) transformation , (43) yieldswith . 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 systemwith and as in (43) and , the second Lyapunov coefficient. Here if there exist critical values of the parameters, wherethen the linear local transformation transforms (47) intowhere . 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,whereThe 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 becomeswhere application of the linear time scaling (41) and the nonlinear time reparametrizationresults in the locally orbitally equivalent systemwith , , 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) aswith . 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 relationswhere the derivatives of with respect to the parameters are given by (49) andThe 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 inBy 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 , whereAs 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 aswe 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 systemwith , 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 equationThen 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 isWe 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.
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.
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.
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.