## Abstract

We discuss the nature of Neimark–Sacker bifurcations occurring in a non-standard numerical scheme, for a class of positivity-preserving systems of ordinary differential equations (ODEs) which undergoes a corresponding Hopf bifurcation. Extending previous work (Alexander & Moghadas 2005*a* *Electron. J. Diff. Eqn. Conf*. **12**, 9–19), it is shown that the type of Neimark–Sacker bifurcation (supercritical or subcritical) may be affected by the integration time-step . The general form of the scheme in the vicinity of a fixed point is given, from which the expression for the first Lyapunov coefficient for two-dimensional systems, valid for arbitrary time-step, is explicitly derived. The analysis shows that this coefficient undergoes an shift with respect to the corresponding coefficient of the original ODE. This could lead to a type of bifurcation which differs from the corresponding Hopf bifurcation in the ODE, due to changes in the sign of the first Lyapunov coefficient as varies. This is especially problematic in the vicinity of certain types of degenerate Hopf bifurcation, at which this coefficient may vanish. We also present a general method to eliminate the possible shift in the bifurcation parameter of the scheme; however, the first Lyapunov coefficient may still be subjected to an shift, leading to a possibly erroneous type of bifurcation. Examples are given to illustrate the theoretical results of the paper with applications to mathematical biology.

## 1. Introduction

The study of the quantitative behaviour of dynamical systems is fundamental to understanding their qualitative dynamics. These systems often depend on parameters that play crucial roles in their qualitative behaviour, in particular with regard to exchange of stability and bifurcations. These, in turn, may influence the observable quantitative behaviour which is reflected in the discrete dynamics of a numerical scheme (NS) that is designed to approximate the solutions of the underlying continuous system (CS). The discrete dynamics will, in general, be influenced by parameters which characterize the NS (e.g. the integration time-step) and which have no counterpart in the CS. These scheme-dependent parameters may well play a significant role in determining changes in discrete dynamics, and should therefore be treated on the same footing as the parameters of the CS. Hence, it is necessary to explore the relationship between the dynamics of the CS and those of the corresponding NS (Brezzi *et al*. 1984; Hofbauer & Iooss 1984; Stuart & Humphries 1996). This relationship is schematically depicted in figure 1, in which CS and NS denote, respectively, the continuous system (usually a set of differential equations) and the discrete system (numerical scheme).

Numerical schemes have largely been employed for the illustration of qualitative dynamics of the CS, frequently through numerical simulations. However, an important question is whether, and to what extent, these simulations faithfully reproduce the dynamics of the underlying CS. Early work in the 1980s (Brezzi *et al*. 1984; Hofbauer & Iooss 1984) and later (Hairer 1994; Hairer & Lubich 1997; Hairer *et al*. 2002) considered the dynamical features, in particular bifurcation behaviour, of widely used standard NSs, including Euler, Runge–Kutta, and predictor–corrector methods (Lambert 1973). For instance, using reduction on the centre manifold to Poincaré normal form, Hofbauer & Iooss (1984) showed the appearance of an invariant limit cycle through a Neimark–Sacker bifurcation in Euler's method with constant integration time-step. They extended this result to include general Runge–Kutta methods. Similar results were obtained by Brezzi *et al*. (1984) for the forward Euler method. It was also shown that the limit cycle in the NS converges uniformly to the limit cycle in the original differential equation, as the integration time-step approaches zero. An alternative approach, backward error analysis, has been discussed for bifurcation and periodicity behaviours of Runge–Kutta methods (Hairer 1994; Hairer & Lubich 1997; Hairer *et al*. 2002). This method formally interprets the NS as an exact solution of a perturbed differential equation.

However, the standard methods considered in these studies have been shown, in general, to undergo ‘scheme-failure’ (Moghadas *et al*. 2003, 2004; Alexander & Moghadas 2005*a*; Moghadas & Alexander in press), in the sense of erroneously representing the qualitative dynamics of the underlying CS with regard to transient, asymptotic and, particularly, bifurcation behaviour. Scheme-failure often occurs when the system parameters, integration time-step, or even initial conditions, are varied in the vicinity of a bifurcation point (Alexander & Moghadas 2005*a*). These may cause a shift in the parameter value at which bifurcation occurs and also, in some cases, a shift in the position of the fixed points corresponding to the critical points of the CS. Notwithstanding these difficulties, these methods also may not preserve some major physical properties modelled by the CS, such as positivity invariance discussed below. This is a common property in physical and biological phenomena that should be inherited by a NS.

To overcome these problems, there has recently been a surge of interest in designing non-standard finite-difference schemes (Mickens 1994, 1999; Moghadas *et al*. 2003, 2004; Dumont & Lubuma 2005; Alexander & Moghadas 2005*a*; Moghadas & Alexander in press), which began to appear in the late 1980s (Patidar 2005). Mickens (1994) provided a formal basis for developing these schemes that were previously applied only to specific problems. The fundamental principle for their construction is to be consistent with the dynamics of the CS (Mickens 2005), in particular to eliminate the shift in positions of the fixed points, and also preserve the positivity invariance of the CS (Alexander & Moghadas 2005*a*). Examples of CSs with the positivity-preserving property include epidemic models in the study of population dynamics (Anderson & May 1991; Alexander & Moghadas 2004, 2005*a*), predator–prey interactions in ecological modelling (Freedman 1980; Moghadas & Alexander in press) and chemical reaction models, such as the Brusselator system (Tyson 1973; Fogler 1999).

Although the bifurcation behaviour of certain standard numerical methods has been investigated for sufficiently small integration step-sizes (Brezzi *et al*. 1984; Hofbauer & Iooss 1984; Hairer 1994; Hairer & Lubich 1997), to the best of our knowledge nothing comparable is available for non-standard NSs. In particular, whether these non-standard schemes are able to capture the actual dynamics of the underlying CSs regardless of the size of integration time-step is an outstanding question, which has not been fully addressed. This is a rather time-consuming and difficult task for general dynamical systems. Therefore, it is appropriate to consider a class of dynamical systems with certain properties, which includes families of models such as those mentioned above. In a series of studies, commencing with Moghadas *et al*. (2004), we introduced the following class of positivity-preserving systems of ordinary differential equations (Alexander & Moghadas 2005*a*):(1.1)where are non-negative real-valued functions. Many models of population dynamics (Freedman 1980; Anderson & May 1991; Alexander & Moghadas 2004, 2005*b*) can be written in the form of (1.1), and we discuss this further in §6.

To study the quantitative dynamics of (1.1), we constructed a non-standard finite-difference scheme of the following general form (Alexander & Moghadas 2005*a*):(1.2)(1.3)for , where is the integration time-step, and and denote the values of at the advanced and current times, respectively. It can be easily shown (by induction) that iterations of the scheme remain in as long as the initial value is chosen in . Thus, the scheme preserves the positivity property of (1.1). Assuming that is a fixed point of (1.2) and (1.3), i.e. at which , and renaming variables by defining as local coordinates centred at , the scheme can alternatively be written as(1.4)It follows immediately from (1.4) that the critical points of (1.1) coincide with the fixed points of the scheme. From now on, we refer to (1.1) as CS and (1.4) as NS.

The relationship between the stability of the fixed points of NS and that of the critical points of CS has been discussed in Alexander & Moghadas (2005*a*). We also established a connection between the dynamics of both CS and NS for codimension-one bifurcations. We showed that when CS undergoes a steady-state bifurcation, in which a simple eigenvalue passes through zero (e.g. pitchfork, transcritical or saddle-node bifurcation), the scheme exhibits a corresponding bifurcation at the same bifurcation parameter value. On the other hand, for a Hopf bifurcation, there is an shift in the value of the corresponding Neimark–Sacker bifurcation parameter of NS. As a consequence, there is also, in general, an shift in the value of the corresponding bifurcation parameter of NS when CS undergoes a Bogdanov–Takens bifurcation. This codimension-two bifurcation appears as a result of the merging of a saddle-node and a Hopf bifurcation (Kuznetsov 1995). We detailed these results for a generalized predator–prey system of Gause-type (Moghadas & Alexander in press) and the Brusselator system (Alexander & Moghadas 2005*a*). However, there are fundamental problems that remain to be addressed.

The aim of this study is to further explore the relationship between the bifurcation behaviour of CS and that of NS, with particular emphasis on Hopf bifurcations. A major concern is to determine whether the nature of a Hopf (Neimark–Sacker) bifurcation in NS would be the same as the corresponding bifurcation occurring in CS. This can be studied by comparing the expressions for the first Lyapunov coefficients (Kuznetsov 1995) of CS and NS. The expressions derived in this paper are valid for arbitrary integration time-step. The problem may become more complicated, however, when CS undergoes a degenerate Hopf bifurcation, as occurs, for example, when the first Lyapunov coefficient vanishes, and therefore the expressions for the second Lyapunov coefficients of CS and NS need to be compared. This paper will be primarily concerned with two-dimensional systems. Generalization to *m*-dimensional (*m*>2) systems requires projecting the dynamics onto the respective centre manifolds of CS and NS, and applying normal form theory on these sub-manifolds (Kuznetsov 1995).

## 2. Notation

Here, we provide the notation that we shall use throughout this paper. For any , , we defineand

As notation for the *α*-derivative, we defineand similarly for . Using this notation, we may write the Taylor expansion of in (1.4) about as(2.1)where , and(2.2)

## 3. Transformation of NS to explicit form

In this section, we shall transform NS into a general explicit form in a local neighbourhood of a fixed point that is required to further analyse the resulting discrete dynamical system. Note that either form of NS in (1.2) and (1.3) or (1.4) is semi-explicit, whose implementation requires an update of variables in advanced time through a Gauss–Seidel type procedure (Varga 2000). The explicit form of NS will ultimately enable us to derive an expression for the first Lyapunov coefficient at the Neimark–Sacker bifurcation for arbitrary integration step-size, and allow comparison with the corresponding coefficient at the Hopf bifurcation in CS. This comparison is needed to determine whether the nature of the bifurcation in NS is consistent with that in CS.

Let . Since , it follows from the implicit function theorem that in a sufficiently small neighbourhood of a fixed point , may be represented as(3.1)where are unknown coefficients to be determined, and we have used the convention of an implied summation from 1 to *m* over repeated pairs of indices (in this case, *j*). Rewriting (1.4) as , substituting (3.1) into this form, and expanding and in their Taylor series about gives(3.2)where when *k*=*j* and 0 otherwise. Note that in (3.2) appears implicitly through , and should be substituted from the expansion given by (3.1). The unknowns in (3.2) can be determined using the following procedure. Having first computed the coefficients , as explained below, we can express for any order explicitly in terms of the , previously derived coefficients from the setand the known values of and . General expressions for for *m*=2 and are given in appendix A, and used in §6 for the simulations of two examples of two-dimensional systems: a predator–prey model and an epidemic model.

In order to determine , we equate first-order terms in (3.2) to obtain(3.3)Defining and using (2.2), it can be seen that, since is arbitrary, equation (3.3) yields(3.4)where .

Alternatively, the expression for in (3.1) can be written as(3.5)where , with . It then follows from (3.4) that , where .

Noting that as , it follows from equation (3.5) that(3.6)From (1.4), it can be seen that as , whose Taylor expansion is given by (see (2.1))(3.7)Comparison of (3.7) with (3.6) leads to , and therefore(3.8)

## 4. Two-dimensional form of NS near Neimark–Sacker bifurcation

For a general two-dimensional system, from (3.1) we have(4.1)where *A* is given by (3.4), with *m*=2. In this case, a simple calculation yieldsIt follows from (3.4) that(4.2)where . Let and be left and right eigenvectors, respectively, of *A* corresponding to the eigenvalue , such that and . Using the normalization , these vectors may be expressed as(4.3)

Using the expression for *A* in (4.2), it follows that(4.4)where . Since the eigenvalues of *A* are the roots of the characteristic equation , using (4.4) it follows that(4.5)where . Note that the expression under the square root is positive sufficiently close to the Hopf bifurcation point and for sufficiently small. An immediate consequence of (4.4) and (4.5) is that, if is an eigenvalue of *J*, then and , so .

Now, we rewrite (4.1) in complex form. Using (4.3), we can transform the real variables and to the complex variable *z* (Kuznetsov 1995) as(4.6)Noting that , we have(4.7)(4.8)By forming the inner product of with (4.1), and denoting as the value of *z* at the advanced time, we obtain(4.9)Substituting (4.7) and (4.8) into (4.9), and defining(4.10)we derive(4.11)This equation can be written in the standard form (Kuznetsov 1995)(4.12)where the expression for may be determined from (4.11) as follows. First, note that(4.13)where is the total order in . Setting and , the right-hand side of equation (4.13) reduces to(4.14)where the limits and are defined such that the arguments of the binomial factors are within the ranges and , respectively. Hence, and . Comparing the expression in (4.14) with in (4.12), setting , and redefining variables leads to the following expression:(4.15)

We emphasize that the expressions for are exact for all . This is important for determining the Lyapunov coefficients below, where their behaviour for arbitrary is of interest.

## 5. First Lyapunov coefficient

Having computed the coefficients in (4.15), we are now in a position to give the normal form of the Neimark–Sacker bifurcation (Kuznetsov 1995, lemma 4.7). Let *β* denote the bifurcation parameter, such that NS undergoes a Neimark–Sacker bifurcation at *β*=0. If is the value of *θ* at bifurcation (*ρ*=1), then provided , , we can transform (4.11) by successive smooth changes of coordinates(5.1)into the normal form(5.2)for *λ* sufficiently close to , where(5.3)

Then, provided , there is a unique closed invariant curve bifurcating from the fixed point as *β* passes through zero. As can be seen from (4.4) and (4.5), the condition , , is satisfied for sufficiently small , since . We now use the expression for to establish a relationship between the first Lyapunov coefficient for the Hopf bifurcation in CS and the corresponding coefficient for the Neimark–Sacker bifurcation in NS.

Since , it follows from (4.2) and (4.10) that

Then from (3.8), (4.5) and (4.10), it can be seen that

Noting that in (4.15) is and , it follows that is also . Letting , it can be seen that , where is defined as the first Lyapunov coefficient in the normal form for the Hopf bifurcation in CS (Kuznetsov 1995, lemma 3.6), obtained from (4.12) using . This implies that is an shift of evaluated at the Hopf bifurcation. Therefore, we have established the following theorem.

*Let* *and* *be first Lyapunov coefficients for CS and NS, respectively, evaluated at the Hopf and Neimark–Sacker bifurcations. Then,*

According to theorem 3.2 of Alexander & Moghadas (2005*a*), there is, in general, an shift in the value of the bifurcation parameter in order for the Neimark–Sacker bifurcation in NS to reproduce the Hopf bifurcation of CS, regardless of the nature of the bifurcation in CS. However, the shift predicted by theorem 5.1 is an independent phenomenon affecting the *type* of bifurcation. For certain systems (e.g. predator–prey model in Alexander & Moghadas (2005*a*) and Moghadas & Alexander (in press)), the shift in the bifurcation parameter of NS disappears, but, in general, an shift may still appear in the value of . These findings highlight the possibility of erroneously reproducing the type of Hopf bifurcation, if the integration time-step is not sufficiently small. Furthermore, if is sufficiently small (i.e. close to a degenerate Hopf bifurcation), the shift in theorem 5.1 may change the sign of for large enough integration time-step, and hence NS would generate a different type of bifurcation (either supercritical or subcritical) to the one predicted by CS.

Even though NS may reproduce the same type of Hopf bifurcation as occurs in CS, there is still an shift in the radius and period of oscillation in the vicinity of the bifurcation. This can be shown as follows. Letting(5.4)and using eqn (3.2) of Alexander & Moghadas (2005*a*), it can be seen that *ξ* satisfies the quadratic , where(5.5)

(5.6)

Taking into account that and at the bifurcation point of CS, it follows that and . Therefore,(5.7)where is the frequency of the oscillation and is the eigenvalue of *J*, at the Hopf bifurcation. Substituting (5.7) into (5.4) gives(5.8)

Now, if we write *z* in the complex form , then and , where describes the rotation angle for the integration time-step in the vicinity of the Neimark–Sacker bifurcation (i.e. ) and . Thus, the period of the oscillation created by the NS in the vicinity of the Neimark–Sacker bifurcation is given by(5.9)where is the period of oscillation in CS. A similar argument shows that the radius of the closed invariant curve undergoes an shift with respect to the radius of the limit cycle in CS (see also Brezzi *et al*. 1984). This behaviour was observed in simulations for the Brusselator system in Alexander & Moghadas (2005*a*).

According to the analysis of this section, the maximum time-step size for which NS is capable of faithfully reproducing the qualitative dynamics of CS may be limited. Moreover, this limitation becomes more stringent as CS approaches an ordinary or degenerate Hopf bifurcation. To remove such restriction and extend the time-step interval without incurring numerical instabilities, previous studies (Mickens 1994; Dumont & Lubuma 2005 and references therein) suggested the use of a more complex form of the time-step, which may also depend on parameters in CS. This function is known as the ‘denominator function’ and has the general form , where *β* denotes the vector of CS parameters (see figure 1). Using first-order perturbation analysis (Alexander & Moghadas 2005*a*), it can be shown that the shift in the Hopf bifurcation parameter and the first Lyapunov coefficient are not affected by the choice of denominator function (*ψ*); however, it may influence terms. These properties can be demonstrated by replacing with *ψ* in the expressions given in appendix A, and expanding them in powers of . In conclusion, the denominator functions are able to eliminate neither the shifts nor the restriction on the maximum time-step in the bifurcation behaviour of the non-standard scheme developed here.

We may take advantage of the structure of the CS equations to revise NS, in an attempt to eliminate shifts in bifurcation behaviour of NS. For example, defining(5.10)(5.11)where is a differentiable, non-negative, real-valued function on , it can be seen that in (1.1) remains unchanged, regardless of the functional form of , for at the Hopf bifurcation. Note that, in general, and will depend on the model parameters *β*, but we suppress this notation here. For a two-dimensional system, by substituting and into (5.6), it can be shown that the coefficient of vanishes, provided(5.12)This provides a degree of freedom in the choice of and to satisfy the condition for the Hopf bifurcation simultaneously in CS and NS, by elimination of the shift in the bifurcation parameter of NS. For an *m*-dimensional system , this shift can be removed through a similar analysis, by applying the first-order perturbation method to equation (3.1) in Alexander & Moghadas (2005*a*). More precisely, by substituting and into eqn (3.13) of Alexander & Moghadas (2005*a*), it follows that the coefficient of in the expansion for the perturbed eigenvalue (in eqn (3.5) of the cited paper) vanishes, provided can be found to satisfy(5.13)where and are, respectively, the left and right eigenvectors of the Jacobian of CS at the Hopf bifurcation point. Note, however, that the bifurcation parameter may still be affected by terms. Elimination of these higher-order terms would require further, and likely more complicated, conditions on .

Remark 5.3 introduces a possible method for revising NS in order to preserve the value of the Hopf bifurcation parameter in CS. However, this method does not, in general, eliminate the shift in the first Lyapunov coefficient, so that theorem 5.1 still holds. This is illustrated in §6*b*. In order to remove this shift, in addition to (5.13), would have to satisfy the algebraic equation resulting fromwhere is defined in §5.

## 6. Applications

In this section, we shall illustrate our results by deriving the first Lyapunov coefficient of NS for two models of real-life dynamical systems in ecological and epidemiological modelling. For the simulations associated with each model, we have used the expressions for and at the bifurcation points, which are given in appendix A.

### (a) Example 1: a predator–prey model

Here, we consider a predator–prey model of Gause-type given by (Freedman 1980)(6.1)where *x* and *y* are the prey and the predator population sizes, respectively; *r*, *k*, *μ* and *D* are positive parameters representing prey intrinsic growth rate, carrying capacity, conversion rate of prey to predator and predator death rate, respectively. The function is called the functional response of predator to prey and satisfies some general assumptions (Moghadas & Alexander in press). For the derivation of equations (6.1) and biological description of the model, the reader may consult May (1976), Freedman (1980) and Murray (1989).

Let , of Holling type II (Holling 1965), which is also the most commonly used type of functional response (Jeschke *et al*. 2004). This model can be written in the form (1.1), with and , where . When , there is a critical point in the first quadrant, which undergoes a supercritical Hopf bifurcation for a certain value of *k* (Moghadas & Alexander in press). Since at , a simple calculation yields that in (5.6), and therefore no shift occurs in the value of the bifurcation parameter of the corresponding NS (Alexander & Moghadas 2005*a*). However, the value of , and consequently the type of bifurcation, may be affected by . Simulation results, depicted in figure 2, demonstrate that although the type of bifurcation is also preserved by NS, the value of increases slightly as increases from 0 to 0.01.

### (b) Example 2: an epidemic model

Let us now consider a susceptible-infected (SI) epidemic model given by Alexander & Moghadas (2004),(6.2)(6.3)where *S* and *I* represent, respectively, the number of susceptible and infected individuals, is the rate of recruitment of individuals (including newborns and immigrants) classified as susceptible, *p* is the proportion of recruited individuals who are vaccinated and therefore completely protected, *μ* represents the natural death rate, *β* is the probability of infection per contact per unit time, *α* is the recovery rate and for *I*, , represents the functional form of a nonlinear force of infection, which satisfies a few general assumptions. The model is biologically well justified and we refer the reader to Alexander & Moghadas (2004) for a detailed description and analysis. It was shown that the model undergoes supercritical and subcritical Hopf bifurcations for certain values of the bifurcation parameter *ν*. This model can be written as (1.1) with and .

Considering the form with (see §4.2 in Alexander & Moghadas (2004)), the condition in (5.6) at the Neimark–Sacker bifurcation of the corresponding NS can be expressed as(6.4)Using the fixed point equation at the bifurcation point (see eqn (2.2) of Alexander & Moghadas (2004)), and defining , it can be seen that *u* satisfies the following equation:(6.5)where . This equation has two positive roots, provided the discriminant of is non-negative, or equivalently (see also inequality (4.7) in Alexander & Moghadas (2004))(6.6)Thus, using , the fixed points of NS at the Hopf bifurcation are given by eqn (2.2) of Alexander & Moghadas (2004), from which and .

We have run the simulations for the fixed points and the corresponding bifurcation parameters of NS, using parameter values given in Alexander & Moghadas (2004). Unlike the previous example, the value of the parameter *ν* at bifurcation undergoes an shift, as shown by the solid curve in figure 3*a*. For these parameter values, CS undergoes a supercritical Hopf bifurcation. However, as depicted in figure 3*b* (solid curve), the sign of the first Lyapunov coefficient changes when , leading to a change in the type of bifurcation in NS from supercritical to subcritical as increases. It is important to note that, in order to apply (5.3), the change in *ν* with must be taken into account as well as the corresponding change in .

We may follow the procedure given in remark 5.2 to eliminate the shift in the value of *ν*, by choosing and , which satisfy (5.12). The results are shown by the dashed curves in figure 3*a*,*b* for *ν* and , respectively. Clearly, while the shift in *ν* disappears, it still causes a change of sign in and in the type of Hopf bifurcation in NS.

## 7. Discussion

Although NSs are essential tools for displaying the qualitative behaviour of dynamical systems, the occurrence of scheme-failure has been a major difficulty requiring a systematic analysis of their structure. Scheme-failure may arise in different ways, the most important being an improper design that may give rise to the loss of some critical properties of the CS (Mickens 2004). If these properties are not inherent to the scheme, it may still fail to reproduce the actual dynamics of the CS, even with small integration time-steps. This has already been shown for the standard adaptive RK45 method, which fails to preserve the positivity property of the original CS (Moghadas *et al*. 2003). This is an important common property of many dynamical systems, in particular in the study of natural phenomena in population dynamics, chemical kinetics, or ecosystem interactions.

In this paper, we discussed the nature of the Neimark–Sacker bifurcation in a non-standard scheme that intrinsically preserves the positivity property of the CS. A relationship between the corresponding types of bifurcation in CS and NS was established by explicitly deriving the first Lyapunov coefficient of NS that is valid for arbitrary integration time-step , and describing its behaviour as a function of . We detailed the results by simulating two examples of a predator–prey model and an SI epidemic model. In the first example, we observed that although agrees in sign with , its value increases as step-length is increased. Not only does this behaviour occur in the second example, but also undergoes a change in sign at a certain (small) step-length, which leads to a transition from supercritical to subcritical Hopf bifurcation. Thus, while NS is capable of producing a Hopf (Neimark–Sacker) bifurcation, it may fail to generate the correct (supercritical) type of bifurcation that occurs in the continuous epidemic model. This clearly shows the importance of theorem 5.1, and emphasizes that reducing the integration time-step to better approximate the value of the bifurcation parameter (as explained in theorem 3.4 of Alexander & Moghadas (2005*a*)) may not always be sufficient to also capture the correct qualitative dynamics of CS. As our analysis has shown, this is particularly problematic in the vicinity of a degenerate bifurcation, for example when is close to zero. This certainly warrants further investigation, by extending the normal form (5.2) to higher order to determine the second Lyapunov coefficient of NS and its behaviour as step-length is varied.

We noted that the expressions for and in CS, and therefore in NS, may not be unique, and we have described a technique using this property for eliminating the shift in the value of the Hopf/Neimark–Sacker bifurcation parameter in *m*-dimensional systems. This non-unique representation of CS and NS affords an opportunity to revise the NS presented here, and whether this can be exploited to provide practical techniques for also eliminating the shift in Lyapunov coefficients (which are especially problematic near degenerate Hopf bifurcations) is a subject of future investigation.

In the present work, we have expressed the general form of NS in the vicinity of a fixed point for *m*-dimensional systems (see (3.1) and (3.2)). This enabled us to derive the general form of NS in the vicinity of a Neimark–Sacker bifurcation for *m*=2 in (4.12). In order to generalize to arbitrary dimensions, we may apply lemma 5.2 of Kuznetsov (1995) to (3.1) and reduce the *m*-dimensional NS to its dynamics on the two-dimensional centre manifold. The analysis can then proceed as discussed in §§4 and 5.

## Acknowledgments

This work was supported in part by the Natural Sciences and Engineering Research Council of Canada (NSERC) and Mathematics of Information Technology and Complex Systems (MITACS). The authors wish to thank the referees for comments that have improved the paper.

## Footnotes

- Received January 16, 2006.
- Accepted March 28, 2006.

- © 2006 The Royal Society