Analytical investigation of single and double Neimark-Sacker bifurcations

The analytical investigation of bifurcations is a very challenging task for many applied scientists and engineers. Often, numerical simulations cannot clarify the complicated dynamics of mechanical systems, in this cases, preprogrammed softwares can be of valid help during the investigation. Also, in the literature, methodology to study bifurcations are presented for most of the cases. However, the presented procedures, are often very hard to be understood from applied scientists with low mathematical background. In this paper we present in details the typical procedure to analyze single and double Neimark-Sacker bifurcations. Especially regarding the double Neimark-Sacker bifurcations of maps, very few sources can be found in the literature, although this kind of bifurcation is very common in many dynamical systems.


Introduction
In the last decades, several authors [2,4,5,9,11] presented methodologies and practical procedures to analytically investigate bifurcations.The typical approach, to analyze a Neimark-Sacker (NS) bifurcation, consists in reducing the order of the system through a center manifold reduction, and then to transform the system into its normal form, in order to analyze the type of motions occurring.This procedure is well known and it has been implemented in several softwares for automatic computation of the bifurcation structure.Nevertheless, in the procedure presented in the literature, many passages are hidden and can be hardly understood from a reader with a weak mathematical background.Especially regarding the double Neimark-Sacker bifurcation, up to our knowledge, the only analytical procedure shown in the literature is in [6], where, although the analysis is considering most of the existing features of a double NS bifurcation and proofs are provided, many basic passages are not explicit and the procedure can be very hard to be understood by engineers or applied scientists.
The aim of this paper is to provide a practical guide for investigating single and double NS bifurcations, according to the most typical approach used in the literature.In the first part of the paper, we present the procedure to investigate a single NS bifurcation, highlighting the passages between the different steps of the analysis that are: transformation to Jordan normal form, center manifold reduction, elimination of nonlinearities through a near identity transformation and transformation to an amplitude map.An analysis of the amplitude map gives information about the behavior of the original system.
In the second part of the paper, we present a procedure to investigate a double NS bifurcation.The main steps are analogous to those of the first part, but the increased dimension of the system makes the procedure much more lengthy.Especially the analysis of the normal form, in case of a double NS bifurcation, can be very complex and very hard to be generalized.

Neimark-Sacker bifurcation 2.1 Mathematical model
We consider a generic map where f = [ f 1 (x; p) ... f n (x; p)] T , x = [x 1 ... x n ] T , p is a scalar real number and n ≥ 2. The trivial solution x 0 = [0 ... 0] T satisfies the equation We consider the case for which the trivial solution x 0 of (1) is stable for p < p cr , while it is unstable for p > p cr .We call p the control parameter of the bifurcation occurring at p = p cr .Assuming that f 1 ..., f n are sufficiently smooth, we expand them in their Taylor series around 0, with respect to x 1 , ..., x n up to the third order, so we can rewrite Eq. ( 1) as where the vector b(x) contains all the nonlinear terms.The stability of the trivial solution depends on the eigenvalues of A(p): the solution is asymptotically stable if and only if all the eigenvalues of A are inside the unite circle of the complex plane, i.e. |µ i | < 1 for i = 1, ..., n.We consider that If the conditions in Eq. ( 4) are satisfied, according to Floquet theory, a Neimark-Sacker bifurcation is occurring for p = p cr [2].

Jordan normal form
Following the procedure indicated in basic textbooks, it is possible to reduce the system to its Jordan normal form, i.e. to reorganize the linear part of the system, in order to decouple the variables related to the bifurcation from the others.
We call H = A| p=p cr and s i , i = 1, ..., n the eigenvectors related to the eigenvalues µ i of H.In the case the eigenvalues µ 3,...,n are real and have algebraic multiplicity equal to 1, we can define the transformation matrix It can be easily verified that Applying the transformation we can rewrite the map in Eq. ( 3) in Jordan normal form If the eigenvalues µ 3−n are not real or have algebraic multiplicity larger than 1, the procedure to obtain the Jordan normal form is slightly different, we refer to [8] for these cases.The matrix T −1 HT controlling the linear part will still be a block diagonal matrix.
In Eq. ( 8), the two variables related to the bifurcation are y 1 and y 2 , and are linearly decoupled from the other variables.

Center manifold reduction
Following the procedure outlined in [2] and [11], we want to reduce the dynamics of the system to its center manifold, i.e. to the two variables related to the bifurcation, y 1 and y 2 .Of course, if the system is already two dimensional, this passage can be skipped.
We define the local center manifold in the form . . .
where t satisfies Eq. ( 8) only for small values of y 1 and y 2 .The cubic terms are neglected, since, after the transformation, they would produce terms higher than the third order.In order to define the coefficients g ihk , we substitute the n − 2 equations of (9) into the first two equations of (8).Then, we substitute these two new equations and the equations in (9) into the remaining n − 2 equations of (8).Collecting terms with the same power order, we obtain 3(n − 2) equations in the 3(n − 2) unknowns g ihk .These equations are organized in a linear system that can be solved in closed form.If now we substitute again the equations in (9) (where g ihk are now known) into the first two equations of (8), we obtain that is the system in Eq. ( 8), limited to its center manifold, i.e. a two dimensional surface in the n dimensional space.Terms higher than the third order, are generated during the transformation and can be neglected.The dynamics of the system in Eq. (10), for small values of y 1 and y 2 , is the same of the system in Eq. (8).

Elimination of nonlinear terms
In order to transform the system into its normal form, we follow the steps outlined in [11] and [4].First of all, we rewrite Eq. (10) in complex form where and ν, α i j ∈ C. Substituting the variables y 1 and y 2 , as expressed in Eq. ( 12), into Eq.( 10), the coefficients ν and α i j can be easily defined as follows ν = µ 2 (13) The next step consists of eliminating all the nonlinear terms not related to internal resonances.This can be done through the near identity transformation.In order to eliminate the second order terms, we apply the transformation to Eq. (11), where hence we obtain where h.o.t.indicates terms higher than the third order.Since we want to eliminate the second order terms, we impose that where the αi j are the coefficient of the third order terms, after the effect of the transformation in Eq. ( 21).So we obtain Collecting terms with the same power order, we obtain so, in order to eliminate the second order terms, we must impose , Then, collecting the coefficients of the third order terms, we obtain Analytical investigation of single and double Neimark-Sacker bifurcations 15 we remind that μ2 = µ 1 and µ 2 µ 1 = 1, so the previous equations can be slightly simplified.The transformation in Eq. ( 21) modifies also the higher order terms, that can be neglected in this procedure.After the transformation, Eq. ( 11) will become With a similar procedure, we can eliminate most of the third order terms.We apply the transformation to Eq. ( 31), where Neglecting terms higher than the third order, we obtain Collecting terms with the same power order, we have In order to eliminate the third order terms, we must impose while, to eliminate the term related to w 2 j w j , we should have ) that has no mathematical sense, since µ 2 − µ 2 2 μ2 = 0, so c 21 → ∞.This is due to the internal resonance between the terms µ 2 w j and α21 w 2 j w j .So we let c 21 = 0.It is possible to generalize the value of the coefficients used for the near identity transformation.While eliminating the second order terms, we had where in case of the third order terms, α hk is substituted by αhk .
After the transformation in Eq. (32), Eq. (31) becomes For the sake of simplicity, from now on we substitute the notation j+1 = f ( j ) with → f ( ).  so, in the vicinity of the bifurcation, Eq. ( 37) can be approximated with the map where e iφ = µ 2 .In order to simplify the calculation, k can be linearized in the following way The next step consists of reducing Eq. ( 39) to an amplitude map.
To do so, we introduce the polar coordinates (r, ψ) ∈ R, where which gives us the following Selecting the absolute value of the map, considering that r ≥ 0, we obtain 2  (44) then, expanding in Taylor series the square root, we have where Instead, selecting the phase of Eq. (43), we have and expanding the arctangent in its Taylor series we have Eq. ( 45) is the normal form of the NS bifurcation.

Bifurcation diagram
The analysis of the bifurcation, is reduced to the analysis of the amplitude map in Eq. ( 45).The trivial solution of Eq. ( 45), corresponds to the trivial solution of Eq. ( 1), while nontrivial solutions of Eq. ( 45), correspond to periodic solutions of Eq. (1).We remind that k > 0 for p > p cr , while k < 0 for p < p cr .
Analyzing Eq. ( 45), it is clear that the trivial solution exists for each value of k and ρ, while it is stable only for k < 0, i.e. for p < p cr .At the same time, to have nontrivial solution, we must have Nontrivial solutions of Eq. ( 45) exist, if and only if k/ρ < 0. So there are two different possibilities: The stability of the nontrivial solution can be analyzed considering the first derivate of Eq. (45 so, the solution is stable for k > 0 and unstable for k < 0, as expected according to standard bifurcation theory.As a practical consequence, we would like to point out that subcritical bifurcations limit the basin of attraction of the trivial solution in the stable region, compromising its robustness and causing unexpected motions, if not properly analyzed.For a direct use of the steps just outlined, in order to define the type of bifurcation occurring, after the center manifold reduction, it is enough to collect the coefficients a hk and b hk in Eq. ( 10) and apply Eqs. ( 14)-( 20), ( 28) and (46), that allow to directly derive ρ.Case studies for this kind of bifurcation are presented in [10] and [3].

Double Neimark-Sacker bifurcation
A double Neimark-Sacker bifurcation is a codimension-2 bifurcation of a fixed point of a discrete time dynamical system, i.e. a fixed point of a map.It occurs in a two dimensional parameter space, when two branches of NS bifurcations are intersecting, as shown in Fig. 2. In correspondence of a double NS bifurcation, two pairs of complex conjugate eigenvalues are on the unit circle of the complex plain, while the other eigenvalues are inside the unit circle.As Fig. 2 shows, in order to have a double NS bifurcation to occur, two parameters should be tuned to the critical value at the same time, this event has zero probability to occur, so our analysis is concentrated in studying what happens in the vicinity of the bifurcation point, where the two NS bifurcations are interacting with each other.Being a finite region of space, it is possible to set the two parameters to be in this region, also in real applications.

Mathematical model
Similarly to the previous section, we consider the generic map where We consider that the stability properties of the trivial solution of the system, are analogous to those shown in Fig. 2. p 1 and p 2 are the control parameters of the bifurcation under study.
The first steps of our analysis are analogous to the ones referred to the single NS bifurcation, we repeat them in this section in order to let the procedure be more understandable.
Assuming that f 1 ..., f n are sufficiently smooth, we expand them in their Taylor series around 0, with respect to x 1 , ..., x n up to the third order, so we can rewrite Eq. (51) as where the vector b 2 (x) contains all the nonlinear terms.In some cases, it may be needed to expand the Taylor series up to the fifth order and keep, during all the procedure, terms up to the fifth order, as explained in [6].We will come back later to this point, during the analysis of the normal form.
The stability of the trivial solution depends on the eigenvalues of A 2 (p 1 , p 2 ): the solution is stable and hyperbolic if and only if all the eigenvalues of A 2 are inside the unite circle of the complex plane, i.e. |µ i | < 1 for i = 1, ..., n, otherwise it is unstable or nonhyperbolic.We consider that for If the conditions in (54) are satisfied, a double NS bifurcation is occurring for (p 1 , p 2 ) = (p 1cr , p 2cr ) [4].

Jordan normal form
As we did in the previous section, we reduce the system to its Jordan normal form, in order to linearly decouple the part of the system related to the bifurcation from the rest of the system.
We call H 2 = A 2 | (p 1 ,p 2 )=(p 1cr ,p 2cr ) and s i , i = 1, ..., n the eigenvectors related to the eigenvalues µ i of H 2 .In the case the eigenvalues µ 5−n are real and have algebraic multiplicity equal to 1, we can define the transformation matrix It can be easily verified that where Applying the transformation we can rewrite the map in Eq. (53) in Jordan normal form As in the previous section, if the eigenvalues µ 5−n are not real or have algebraic multiplicity larger than 1, the procedure to obtain the Jordan normal form is slightly different (see [8]), but the matrix T −1 2 H 2 T 2 , controlling the linear part, will still be a block diagonal matrix.In Eq. ( 59), the variables related to the bifurcation are (y 1 , y 2 , y 3 , y 4 ), and are linearly decoupled from the other variables.

Center manifold reduction
If the dimension of the system is larger than 4 (n > 4), we can reduce it to 4 through the center manifold reduction, as we showed in the case of a single Neimark-Sacker bifurcation.The procedure to be used in the case of a double NS bifurcation is very similar to the one already shown, with the difference that the center manifold is now a 4 dimensional subspace of the n dimensional space, so it has much more coefficients than in the previous case.The general form of the center manifold, approx-imated to the second order terms, is as follows 1, j + g 50200 y 2 2, j + g 50020 y 2 3, j + g 50002 y 2 4, j +g 51100 y 1, j y 2, j + g 51010 y 1, j y 3, j + g 51001 y 1, j y 4, j +g 50110 y 2, j y 3, j + g 50101 y 2, j y 4, j + g 50011 y 3, j y 4, j . . .g n2000 y 2  1, j + g n0200 y 2 2, j + g n0020 y 2 3, j + g n0002 y 2 4, j +g n1100 y 1, j y 2, j + g n1010 y 1, j y 3, j + g n1001 y 1, j y 4, j +g n0110 y 2, j y 3, j + g n0101 y 2, j y 4, j + g n0011 y 3, j y 4, j In order to define the 10(n−4) coefficients we substitute the n−4 equations of (60) into the first four equations of (59).Then, we substitute these four new equations and the equations in (60) into the remaining n − 4 equations of (59).Collecting terms with the same power order, we obtain 10(n − 4) equations in the 10(n − 4) unknowns g ihklm .These equations are organized in a linear system that can be solved in closed form.Substituting again the equations in (60) into the first four equations of (59), we obtain h+k+l+m=2,3 a hklm y h 1, j y k 2, j y l 3, j y m 4, j h+k+l+m=2,3 b hklm y h 1, j y k 2, j y l 3, j y m 4, j h+k+l+m=2,3 c hklm y h 1, j y k 2, j y l 3, j y m 4, j h+k+l+m=2,3 d hklm y h 1, j y k 2, j y l 3, j y m that is the system under study, limited to its center manifold.The dynamics of the system in Eq. ( 61) is the same of the system in Eq. ( 51), for small values of (y 1 , y 2 , y 3 , y 4 ).As we did in the previous section, from now on we substitute the notation j+1 = f ( j ) with → f ( ).

Elimination of nonlinear terms
Similarly to the case of a single NS bifurcation, we rewrite the system in complex form, according to the change of variables the system in Eq. (61) becomes where α hklm , β hklm ∈ C. Substituting the variables (y 1 , y 2 , y 3 , y 4 ), as expressed in Eq. ( 62), into Eq.( 61) we can define the values of the coefficients α hklm and β hklm .For the second order terms we have the coefficients Regarding the coefficients of the third order terms, we write here only those that will not be eliminated in the next passages, i.e. α 2100 , α 1011 , β 0021 and β 1110 .These are The next step consist in eliminating all the nonlinear terms not related with internal resonances.We apply the following near identity transformation to Eq. ( 63), where As in the case of a single NS bifurcation, choosing properly the values of the coefficients e hklm and f hklm , all the second order terms can be eliminated.Of course, this procedure will modify the terms higher than the second order.The procedure is similar to the one shown in the previous section, but in this case it is much more lengthy due to the higher dimension of the system.For this reason, we skip this passage and we write directly the values of the coefficients, which follow the same rule of Eq. (36), i.e.
As a result of this transformation we will obtain the system As we did before, we write only the values of the coefficients related to the bifurcation, that are With an analogous transformation, we eliminate the third order terms not related to internal resonances, i.e. we apply the transformation to Eq. ( 82), where The values of the coefficients in the transformation are This procedure allow us to eliminate all the third order terms, except the ones related to w 2 1 w1 and w 1 w 2 w2 in the first equation of (82), and related to w 2 2 w2 and w 1 w1 w 2 in the second equation of (82).These terms are related to internal resonances and cannot be eliminated.From Eqs. ( 89)-( 90) we can easily identify the terms related with internal resonances, in fact, to eliminate those terms, the coefficient e hklm should tend to infinit.For example, to eliminate the term related to w 1 w 2 w2 , we should have After this lengthy procedure, and considering that w w = |w| 2 , we obtain the following system

Reduction to an amplitude map
We introduce two parameters so, in the vicinity of the bifurcation, Eq. ( 92) can be approximated with where e iφ 2 = µ 2 and e iφ 4 = µ 4 .In order to simplify the following calculation and have a direct relationship between k 1 , k 2 and the bifurcation parameters p 1 , p 2 , we linearize k 1 , k 2 in the following way It is possible to reduce the map in Eq. ( 94) to an amplitude map.
To do so, we introduce the polar coordinates (r 1 , r 2 , ψ 1 , ψ 2 ) ∈ R, where w 1 = r 1 e iψ 1 , w 2 = r 2 e iψ 2 . (96) Applying the transformation in Eq. (96) to Eq. ( 94), we have Separating the absolute value from the phase, considering that r 1 , r 2 > 0, we have expanding the square roots in Taylor series we obtain where We now consider the phase of the map in Eq. ( 97), so we have expanding the arctangent in its Taylor series we obtain where

Analysis of the normal form
It is possible to investigate the effect of the double NS bifurcation, analyzing the fixed points of Eq. (99).If there are fixed points different from the trivial solution, the map in Eq. (51) admits periodic solutions.If the fixed point is semitrivial, i.e. it lays on one axis, the map in Eq. (51) will show a periodic motion, while, if the fixed point is nontrivial, the map in Eq. (51) will show a quasiperiodic motion.According to the scheme: fixed point (r 1 , 0) → periodic motion with amplitude r1 and frequency depending on Eq. ( 103 The stability properties of these motions are analogous to the stability properties of the fixed points of Eq. ( 99).
There exists a fixed point on the r 1 axis (r 1 , 0) for so, in order to have a fixed point on the r 1 axis, we must have k 1 /a 11 < 0. The semitrivial fixed point is stable if and only if the eigenvalues of the Jacobian matrix have absolute value less than 1, i.e.
If only one of the two inequalities is verified, the fixed point is a saddle, if both are not verified it is totally unstable.Considering that r1 = √ −k 1 /a 11 and the conditions for stability in Eq. ( 110), if a 11 > 0, this semitrivial fixed point is necessarily unstable, otherwise its stability depends also on a 21 and k 1 .
We now verify the existence of semitrivial fixed points on the r 2 axis (0, r2 ).There exists a semitrivial fixed point on the r 2 axis for similarly to the previous case, in order to have the semitrivial fixed point, we must have k 2 /a 22 < 0. The Jacobian matrix in correspondence of (0, r2 ) is If only one of the two inequalities is verified, the fixed point is a saddle, if both are not verified it is totally unstable.As in the previous case, considering that r2 = √ −k 2 /a 22 and the conditions for stability in Eq. (??), if a 22 > 0, this semitrivial fixed point is necessarily unstable, otherwise its stability depends also on a 12 and k 2 .
There exists a general fixed point (r In order to have real solution we must have r * 2 1 > 0 and r * 2 2 > 0. To analyze the stability of this fixed point, we have to study the eigenvalues of the Jacobian matrix If both the eigenvalues are inside the unit circle of the complex plane, the nontrivial fixed point is stable, and it corresponds to a stable quasiperiodic motion of the map in Eq. ( 51).If only one eigenvalue is out of the unit circle, the nontrivial fixed point is a saddle, so the corresponding quasiperiodic motion of the map in Eq. (51) is unstable.If both the eigenvalues are out of the unit circle, the fixed point is a repellor.As better explained in [6], the nontrivial solution (r * 1 , r * 2 ) arises from a pitchfork bifurcation of one of the two semitrivial solutions (r 1 , 0) or (0, r2 ).
In Fig. 3, we show an example of a possible bifurcations structure in the vicinity of a double NS bifurcation.In the figure, the two single NS bifurcations are both supercritical.In regions B and F there is no interaction between the two bifurcations.In