Stability analysis of degenerate Hopf bifurcations for discrete-time systems

M. B. D'Amico¹, E. E. Paolini and J. L. Moiola^1,2

Departamento de Ingeniería Eléctrica y de Computadoras, Universidad Nacional del Sur
Avda. Alem 1253 - B8000CPB Bahía Blanca, Argentina
mbdamico@criba.edu.ar
¹ CONICET (Consejo Nacional de Investigaciones Científicas y Técnicas)
² Mathematical Institute, University of Cologne, 50931 Cologne, Germany

Abstract ¾ A methodology for the stability analysis of invariant cycles emerging from Hopf bifurcations in discrete-time nonlinear systems is presented. The technique is formulated in the so-called frequency-domain and it is based on the Nyquist stability criterion and a higher-order harmonic balance method. The study of a planar cubic map is included for illustration.

Keywords ¾ Hopf Bifurcation. Discrete-Time Nonlinear Systems. Frequency-Domain. Harmonic Balance Method. Stability Index.

I. INTRODUCTION

The Hopf bifurcation theorem (HBT) for maps describes the appearance of an invariant cycle when one parameter of the system is varied appropriately. Assuming that the fixed point changes its stability, the emerging bifurcation can be supercritical or subcritical denoting the birth of stable or unstable cycles for parameter values larger or smaller than the critical one, respectively. This behavior is similar to that observed in continuous-time nonlinear systems as well as in time-delayed nonlinear systems. Consequently, a technique formulated in the frequency-domain for single-input single-output (SISO) and multiple-input multiple-output (MIMO) discrete-time systems has been introduced in D'Amico et al. (2002) to deal with this characteristic bifurcation. The formulas capture the dynamical behavior of the emerging invariant cycle using concepts from control theory and a second-order harmonic balance method. These results are extensions of the earlier developments obtained by Allwright (1977), and Mees and Chua (1979) for continuous-time systems.

To have a better approximation of the cycle or the possibility of studying more complex dynamical structures, it is necessary to use higher-order expansions of the classical Hopf normal form. This extension is analogous to consider the higher-order approximations obtained using the harmonic balance method in the frequency-domain. However, some care should be exercised when this result is applied to the discrete-time case (Robinson, 1999) as the nonlinear maps frequently exhibit additional dynamical phenomena, such as weak and/or strong resonances. We will not address this issue on this paper, and we will concentrate on deriving the higher-order approximation of the emerging invariant cycle, and on developing algebraic expressions of the so-called stability indices to establish the stability of the cycle even in degenerate Hopf bifurcations (Iooss, 1979; Shilnikov et al., 2001). These indices allows the comparison of the results obtained via the frequency-domain approach with those given by the classical normal form method (Whitley, 1983; Glendinning, 1994; Balibrea and Valverde, 1999). The approximation of the invariant cycle is based on a higher-order harmonic balance resembling the procedures followed by Mees (1981) and Moiola and Chen (1996) for continuous-time systems.

The conditions for detecting degenerate Hopf bifurcations may be translated to the discrete-time case using the frequency-domain approach. Moreover, some of the results can be applied to a much more complex theoretical construction, such as the Poincaré map (Kuznetsov, 1995), to study the stability of quasiperiodic motion in continuous-time nonlinear circuits and systems (see, for instance, Bi and Yu, 1999).

The paper is organized as follows. In Section II, higher-order formulas to determine the stability of the invariant cycle emerging from a Hopf bifurcation are derived. The study of a planar cubic map near a degenerate condition is presented in Section III. Finally, in Section IV some conclusions are given.

II. HOPF BIFURCATION IN THE FREQUENCY-DOMAIN

Let us consider the discrete-time nonlinear system

(1)

where x_k, , , , , is the iteration index, is the bifurcation parameter and is a smooth (C^2q+1, q ³ 1) function. All the matrices may have explicit dependence on m, and A may be the zero matrix.

Many distinct but equivalent feedback representations for (1) can be obtained by introducing an arbitrary matrix (which may also depend on m) and then applying the z-transform. Thus,

	(2)
	(3)

This representation suggests that (1) can be seen as a feedback interconnection between the linear transfer matrix and the nonlinear function , defined by (2) and (3), respectively. The equivalent feedback system is presented in Fig. 1, where d_k represents perturbations, noise effects, etc., v_k is an external reference input, u_k is the control variable and y_k is the output. Observe that d_k and v_k are set to zero as in the continuous-time version of the HBT.

Local dynamical behavior is analyzed by means of the linearization of the open-loop system, given by G(z; m) J(m) where is the Jacobian matrix. As a result, the crossing of a simple pair of complex eigenvalues of (1) through the unit circle for a given value m = m₀ is equivalent to the crossing of one eigenvalue of , denoted as , over the critical point -1 + i0 for certain values w₀ and m₀.

A. Higher-order harmonic balance method

A frequency-domain approach to analyze Hopf bifurcations for MIMO discrete-time systems has been introduced in D'Amico et al. (2002). The technique is based on the application of a second-order harmonic balance to capture the periodic solution emerging from the bifurcation. A natural way to obtain a better approximation is the addition of higher-order terms via a higher-order harmonic balance method.

Let us fix m close to m₀, so the Nyquist diagram of lies near the critical point -1 + i0. In this case, if a periodic solution exists, it can be written as

(4)

with and . Expanding the nonlinear function with respect to e_k in Taylor series up to the (2q + 1)-order and replacing e_k with (4), we obtain

where the coefficients F^r depend on the vectors E^r. Assuming also that the harmonics E^r are function of the amplitude q of the periodic solution,

(5)

with V_rj as intermediate variables grouping the contributions of E^r and j varying from r (or 2 for r = 0) in steps of 2 up to 2q + 1 (or 2q for r even), the coefficients F^r can be expressed as

(6)

where j varies in the same way as before and each W_rj is a function of the higher-order derivatives evaluated at and the vectors V_r'j' with r' £ r and r' + j' £ r + j except for the r = 0 case (see Mees, 1981 for details).

Figure 1: An equivalent feedback representation of (1).

The harmonic balance equations, obtained equating the input and output of the linear part of the system, are given by

(7)

Substituting (5) and (6) into (7), we have

(8)

To avoid resonances, we will assume that the matrix is nonsingular for r > 1. In that case, (8) can be rewritten in the compact form

with and then each V_rj may be calculated if the vectors V_1j' with j' £ j are known.

For r = 1, is singular so that it is not possible to solve (8) directly. To overcome this difficulty, we suppose V₁₁ = v, where v is the normalized right eigenvector of associated with the eigenvalue , and V_1j, with j = 3, 5, ..., are orthogonal to v. Thus, each vector V_1j is obtained from

where j = 1, 2, ..., q and P = I - vv^T is the projection on the subspace orthogonal to V₁₁. The idea behind these assumptions is that q fixes the amplitude of E¹ in the direction of v, acting the remainder terms as corrections in directions orthogonal to that of the main contribution. Then, the harmonic balance equation for r = 1 is

	(9)

Premultiplying both sides by u^T, which is the normalized left eigenvector of associated with , and assuming that

(10)

then (9) can be expressed as

Finally, equating the terms of the same power in q,

It is clear that all the x_j expressions can be calculated if the vectors V_1,2j+1 are known.

From a control theory viewpoint, Eqn. (10) is interpreted as the intersection between the Nyquist locus of and a curve depending on q starting at -1 + i0. If this intersection occurs at q_R = then q_R and w_R are the amplitude and frequency of the periodic solution, respectively.

B. Algebraic expression of the stability indices

Based on the graphical approach presented previously, an analysis for small perturbations of the point q_R may reveal the stability of invariant cycles in maps (Mees and Chua, 1979). Another alternative is the calculation of the stability indices (or curvature coefficients). The computation of the first-order stability index through the formulation of a second-order harmonic balance has been presented in D'Amico et al. (2001). In this section, we will derive algebraic expressions of higher-order indices using the frequency-domain approach. The interested readers can obtain more details in Moiola and Chen (1996) for the analogous continuous-time case.

Assuming that the matrix G(z; m) J(m) possesses an eigenvalue it is easy to verify that the matrix has the same eigenvalue. Moreover, at the criticality, we have and , which are actually the two corresponding complex eigenvalues of (1) of the discrete-time version of the HBT for a given value m = m₀.

Now, let us suppose that when the parameter m is larger than m₀, the pair of complex eigenvalues of (1) crosses the unit circle taking the value , with r > 1. Furthermore, let us consider that under this condition the analysis in the frequency-domain ensures a periodic solution with an approximate frequency w and small amplitude q. Therefore, we can write

(11)

where G'(z; m) and G''(z; m) are the first and second-order derivatives of G(z; m) with respect to z, respectively.

On the other hand, substituting (5) into (9) for r = 1 and considering that V₁₁ = v, W_1,2j+1 and V_1,2j+1 with j = 1, 2, ..., q are known,

Premultiplying both sides by u^T, replacing with (11) and extending the result in terms of the vectors V_1,2j+1 and p_j,

(12)

where v' = dv/dz and and for j = 1, 2, ..., q.

Then, if q is close to zero, the first-order approximation is given by

(13)

By substituting (13) into (12) and grouping together the coefficients of equal power in q,

with h = u^TG'(z; m)J(m)v, which can be written in the compact form

(14)

Observe that, instead of using (13) as a first-order approximation of , it can be used the more accurate Eqn. (14). In fact, following a similar reasoning with (12) and (14), it is possible to calculate the next higher-order approximation. This procedure can be continued in the same way until obtaining the desired order.

Equation (14) allows us to determine the stability of the invariant cycle emerging from a Hopf bifurcation. Defining and considering the real part of (14), we have

(15)

Taking into account that r > 1, the left-hand side of this equation is always greater than zero. Therefore, to have a solution for small q² > 0, it is necessary that

with p₁ calculated as in Moiola and Chen (1996). Although s₁ depends on (the frequency of the exact periodic solution) a reasonably accurate (local) approximation can be obtained computing s₁ at w₀, i.e. the frequency at which the Nyquist eigenlocus of passes over the critical point -1 + i0.

In the case that s₁ = 0 at criticality, the stability of the emerging invariant cycle can be determined analyzing the coefficient of the term q⁴ in the expansion (15). Defining the second-order stability index as , there will exist a periodic solution only if s₂ > 0. In a similar way, if both s₁ and s₂ vanish for certain critical combination of the system parameters, it will be necessary to study the sign of the coefficient corresponding to the term q⁶, and so on.

III. EXAMPLE

Let us consider the planar cubic map

(16)

where d₁, d₂ are constants which do not vanish simultaneously and a, b are bifurcation parameters. Notice that if d₁ and d₂ were equal to zero, the nonlinear terms in (16) would disappear so that the planar map would become linear, acting as a center.

One of the fixed points of the cubic system is (x¹, x²) = (0, 0) and its eigenvalues in the time-domain are given by For r² = a² + b² < 1 the origin is stable and as r² is increased to 1, the pair of complex conjugated eigenvalues crosses the unit circle from inside to outside changing the stability of the fixed point. In fact, the planar map exhibits a Hopf bifurcation at the critical point r₀ = 1.

In order to determine the stability of the invariant cycle emerging from the criticality r = r₀, (16) can be transformed into (2) and (3) considering

Then,

with m as the parameter vector (a, B). Linearizing f(e_k; m) around the fixed point , we obtain J(m) = b/d₁ and thus,

The only eigenvalue of this SISO feedback system is , and then the right and left eigenvectors are given trivially by u^T = v = 1. Moreover, V₁₃ = 0, and in general for k = 4, 5, ... . Therefore,

As a result, and the first-order stability index at the criticality is

(17)

As can be seen, the sign of s₁, and thus the stability of the invariant cycle, depends on the values of f, d₁ and d₂. This fact can be corroborated via numerical simulations, fixing f = 0.515, d₁ = 2 and using d₂ as a control parameter. Figure 2 shows the invariant cycle obtained for d₂ = 0.7 and r = 1.015 > 1. In this case, the stability index is s₁ = 0.1408 > 0 so that the invariant cycle is stable. On the other hand, for d₂ = 1.4 and r = 0.991 < 1, s₁ = -0.088 and thus the emerging cycle is unstable (Fig. 3). For values of r beyond to r₀, these cycles interact with the other two fixed points of the system, and disappear.

Figure 2: Stable invariant cycle obtained for d₁ = 2, f = 0.515, r = 1.015 and d₂ = 0.7 (s₁ = 0.1408).

The planar map (16) develops a degenerate Hopf bifurcation if the index s₁ vanishes at the criticality, i.e. if

d₂ = d₁ tan f (18)

Therefore, to establish the stability of the emerging invariant cycle the computation of the next higher-order index is required.

For a fourth-order harmonic balance, it is found that V₀₄ = V₂₄ = 0,

and then the second-order stability index at the critical point is

Taking into account (18), this expression can be rewritten as

(19)

Since s₂ > 0 for all values of d₁, the cycles emerging from the degenerate bifurcation (s₁ = 0) are stable. Once again, numerical results verify this prediction. As before, we choose f = 0.515, d₁ = 2 and vary d₂ and r. The invariant cycle obtained for d₂ = 1.1302 and r = 1.0001 is presented in Fig. 4. In this case, d₂ is very close to the value d_2d = 1.131547 corresponding to the degenerate condition and thus, s₁ = 4.399 ´ 10^-4. Although r is closer to the unit circle than in the situation of Fig. 2, the amplitude of the invariant solution is larger. This behavior becomes more noticeable as d₂ is closer to d_2d. For values of d₂ slightly larger than d_2d, we should observe a connection between both stable and unstable cycles that occurs generally through a limit point bifurcation. However, this effect is not appreciated in this system because of the disappearance of the cycles explained previously. This phenomenon of cycle connection seems to be common in discrete-time systems of dimension equal or greater than 2, certainly related to global bifurcations and resonance phenomena, and thus explain the interest to study these behaviors in the specialized literature (Iooss, 1979; Whitley, 1983; Shilnikov et al., 2001).

Figure 3: Unstable invariant cycle obtained for d₁ = 2, f = 0.515, r = 0.991 and d₂ (s₁ = -0.088).

It is worth noticing that the expressions of the first and second stability indices, given by (17) and (19), respectively, are in complete agreement with those obtained applying the normal form technique (Iooss, 1979).

IV. CONCLUSIONS

A frequency-domain approach to detect Hopf bifurcations for discrete-time systems using a higher-order harmonic balance method has been presented, and algebraic expressions that determine the stability of the emerging invariant cycle have been derived. The application of the main results is shown studying a planar cubic map.

It is noticeable that the coalescence between stable and unstable invariant cycles has not been detected in the system, although an extensive numerical search has been performed in the parameter space. However, this fact confirms one of the first remarks in the introduction: the analysis of Hopf bifurcations in discrete-time systems is more difficult than that for continuous-time systems. A relationship between the coalescence of invariant cycles and the effects of resonances will be pursued in the forthcoming investigations following the preliminary results obtained by Frouzakis et al. (1991) and Peckham et al. (1995), but including the information of higher-order bifurcation formulas.

Figure 4: Stable invariant cycle for d₂ = 1.1302 near the degenerate condition d_2d.

Acknowledgments
M. B. D. and J. L. M. appreciate the support of CONICET and a grant from SGCyT of the UNS. J. L. M., in addition, acknowledges the financial support of the Alexander von Humboldt Foundation and the hospitality of the University of Cologne.

REFERENCES
1. Allwright, D. J., "Harmonic balance and the Hopf bifurcation theorem", Math. Proc. of Cambridge Phil. Soc. 82, 453-467 (1977).         [ Links ]
2. Balibrea, F. and J. C. Valverde, "Bifurcations under nondegenerated conditions of higher degree and a new simple proof of the Hopf-Neimark-Sacker Bifurcation Theorem", J. of Math. Analysis and Applications 237, 93-105 (1999).         [ Links ]
3. Bi, Q. and P. Yu, "Double Hopf bifurcations and chaos of a nonlinear vibration system", Nonlinear Dynamics 19, 313-332 (1999).         [ Links ]
4. D'Amico, M. B., J. L. Moiola and E. E. Paolini, "Hopf bifurcation for maps: A frequency domain approach", IEEE Trans. on Circuits and Systems-I 49, 281-288 (2002).         [ Links ]
5. D'Amico, M. B., E. E. Paolini and J. L. Moiola, "On the stability of Hopf bifurcation for maps", Proc. Workshop NDES'2001, Delft, The Netherlands, 209-212 (2001).         [ Links ]
6. Frouzakis, C. E., R. A. Adomaitis, and I. G. Kevrekidis, "Resonance phenomena in an adaptively-controlled system", Int. J. of Bifurcation and Chaos 1, 83-106 (1991).         [ Links ]
7. Glendinning, P., Stalility, Instability and Chaos: An Introduction to the Theory of Nonlinear Differential Equations, Cambridge University Press, New York (1994).         [ Links ]
8. Iooss, G., Bifurcations of Maps and Applications, Mathematical Studies 36, North Holland, Amsterdam (1979).         [ Links ]
9. Kuznetsov, Y. A., Elements of Applied Bifurcation Theory, Springer Verlag, New York (1995).         [ Links ]
10. Mees, A. I., Dynamics of Feedback Systems, John Wiley & Sons, Chichester (1981).         [ Links ]
11. Mees, A. I. and L. O. Chua, "The Hopf bifurcation theorem and its applications to nonlinear oscillations in circuits and systems", IEEE Trans. on Circuits and Systems 26, 235-254 (1979).         [ Links ]
12. Moiola, J. L. and G. Chen, Hopf Bifurcation Analysis: A Frequency Domain Approach, World Sci. Series on Nonlinear Science A 21, World Scientific Publishing, Singapore (1996).         [ Links ]
13. Peckham, B. B., C. E. Frouzakis and I. G. Kevrekidis, "Bananas and banana splits: A parametric degeneracy in the Hopf bifurcation for maps", SIAM J. Math. Anal. 26, 190-217 (1995).         [ Links ]
14. Robinson, C., Dynamical Systems-Stability, Symbolic Dynamics, and Chaos, CRC Press, Boca Raton, USA (1999).         [ Links ]
15. Shilnikov, L. P., A. L. Shilnikov, D. Turaev and L. O. Chua, Methods of Qualitative Theory in Nonlinear Dynamics-Part 2, World Scientific Publishing, Singapore (2001).         [ Links ]
16. Whitley, D., "Discrete dynamical systems in dimensions one and two", Bull. London Math. Soc. 15, 177-217 (1983).         [ Links ]