Prediction of homoclinic bifurcation: the elliptic averaging method

Chaos, Solitons and Fractals 11 (2000) 2251±2258

www.elsevier.nl/locate/chaos

Prediction of homoclinic bifurcation: the elliptic averaging method M. Belhaq *, F. Lakrad Group in Nonlinear Oscillations and Chaos, Faculty of Sciences Aõn Chock, Laboratory of Mechanics, Maarif BP 5366, Casablanca, Morocco Accepted 28 July 1999

Abstract A criterion to predict bifurcation of homoclinic orbits in strongly nonlinear autonomous oscillators is presented. The averaging method combined formally with the Jacobian elliptic functions is applied to determine an approximation of limit cycles near homoclinicity. We then introduce a criterion for predicting homoclinic orbits, based on the collision between the bifurcating limit cycle and the saddle equillibrium. In particular, we show that this criterion leads to the same results as the standard Melnikov technique. Explicit applications of this criterion to quadratic nonlinearities are included. Ó 2000 Elsevier Science Ltd. All rights reserved.

1. Introduction The last forties saw a great development of perturbation methods, based on trigonometric functions, for approximating solutions to weakly nonlinear oscillators in the form x ‡ c1 x ˆ eg…l; x; x_ †:

…1†

Here c1 is a positive constant, e a small positive parameter, g a polynomial function of its arguments and l is a control parameter. Classical methods, such as harmonic balance, Lindstedt±Poincare, Krylov±Bogolioubov±Mitropolski, and multiple scales [1±7], have been conducted to approximate periodic solutions of Eq. (1). Recently, attention was devoted to study solutions of strongly nonlinear oscillators in the form x ‡ c1 x ‡ c2 f …x† ˆ eg…l; x; x_ †:

…2†

Here c1 and c2 are ®xed constants, f …x† includes cubic or quadratic polynomial terms and g…l; x; x_ † is an arbitrary nonlinear function of its arguments. Most of the above classical approach has special attention have been extended by introducing the Jacobian elliptic functions [8±21]. Indeed, there was unanimity in the qualitative improvements of the approximation given by such functions in comparison with the trigonometric ones [8±21]. In particular, the use of Jacbian elliptic functions gives an excellent approximation of the periodic orbits even near the separatrices just prior to connection. For instance, Barkham and Soudack [8,9], Soudack and Barkham [10,11] and Yuste and Bejarano [12] used the Krylov±Bogolioubov method to provide approximate solutions of a strongly nonlinear oscillator in terms of Jacobian elliptic functions. Bejarano et al. [13] and Yuste Bejarno [14] used the elliptic functions to approximate periodic solutions in *

Corresponding author. E-mail address: [email protected] (M. Belhaq).

0960-0779/00/$ - see front matter Ó 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 9 6 0 - 0 7 7 9 ( 9 9 ) 0 0 1 4 4 - 7

2252

M. Belhaq, F. Lakrad / Chaos, Solitons and Fractals 11 (2000) 2251±2258

nonlinear oscillators. Rand [15] applied the averaging method to quadratic f …x†. He observed that the approximations based on elliptic functions give much better results, when compared to the traditional trigonometrical techniques. Yuste and Bejarano [16] extended the methods of harmonic balance to a certain class of nonlinear oscillators by introducing Jacobian elliptic functions. Garcia-Margallo and Bejarano [17] used generalized Fourier series and elliptic functions to determine solutions to ®rst order using harmonic balance. Coppola and Rand [18,19] applied the averaging method to cubic f …x† with c1 and c2 slowly varying and used Macsyma code to implement elliptic functions. For other applications of elliptic functions see [20,21]. In this paper, we combine the averaging method with Jacobian elliptic functions and apply the collision criterion of homoclinicity given in [22,23] to obtain an analytical approximation of homoclinic bifurcations that occurred in Eq. (2). The classical, mathematically rigorous approach to predict such bifurcations is the Andronov±Melnikov method. Geometrically, it routinely approximates vanishing distance between the separatrices. Recently, Belhaq and Fahsi [24] proposed an approach based on trigonometric averaging technique, formally approximating in®nite period of the bifurcating periodic orbits. Another criterion based on the collision, at the homoclinic bifurcation, of the periodic orbit with the saddle equilibirium is presented in [22,23]. In all these approaches [22±24], however, the periodic solutions are approximated using perturbational techniques with trigonometric functions. More recently, Belhaq et al. [25] proved for general strongly nonlinear oscillators that the prediction obtained formally by combining the collision criterion and the Jacobian elliptic Lindstedt±Poincare method gives the same results as the Melnikov method. This collision criterion was successfully applied to predict homoclinic bifurcation in a three-dimensional system [26]. Here the approximation of multiple scales technique was conducted to construct a higher order expansion of the bifurcating periodic solution using trigonometric function. The collision criterion was then applied and a critical parameter of homoclinicity was also obtained. The object of this paper is to take advantage of the results presented in [25] based on the elliptic Lindstedt±Poincare technique to derive a condition of homoclinity using the elliptic averaging method. We show that our results coincide with those given by the Melnikov method to leading order. For illustration, we include explicit calculations for quadratic nonlinearities f …x†. 2. Homoclinic collision criterion Consider the unperturbed generating equation of (2) x ‡ c1 x ‡ c2 f …x† ˆ 0:

…3†

For f …x† ˆ x2 or f …x† ˆ x3 , Eq. (3) has an exact analytical solution in terms of Jacobian elliptic functions. Assume that Eq. (3) has a homoclinic orbit C0 to a saddle point S…as ; bs †. Suppose that for small e > 0, an isolated periodic solution survives in the vicinity of C0 and may bifurcate from a homoclinic connection Ce near C0 for ®xed e at some critical parameter value lc . The classical Andronov±Melnikov method [27,28] to predict such a bifurcation is based on the so-called splitting function b. This function is de®ned by considering a one-dimensional local cross-section R to the stable manifold W s and de®ning a coordinate f along R such that f ˆ 0 corresponds to the point of intersection with W s . The splitting function b ˆ fu denotes the f-value of the intersection of W u with R (see Fig. 1). Therefore, the condition for the homoclinic bifurcation to occur is given by b ˆ 0:

…4†

As a variation of the theme, another criterion proposed here is based on the distance between the bifurcating periodic solution and the hyperbolic saddle point S ˆ S…as ; bs †. To be more speci®c, let x…t† ˆ …x…t†; x_ …t††

…5†

be an approximation of the periodic orbit of the perturbed equation (2), located in the vicinity of the homoclinic oribit C0 . Denote by X …x…A; l†; x_ …A; l†† the coordinates of the intersection point of the periodic

M. Belhaq, F. Lakrad / Chaos, Solitons and Fractals 11 (2000) 2251±2258

2253

Fig. 1. Homoclinic bifurcation according to the Andronov±Melnikov method. Fig. 2. Homoclinic bifurcation according to the collision criterion.

solution and an axis (D), connecting the saddle S to a focus F interior to the periodic orbit (see Fig. 2). In the limit l ! lc , the condition to be satis®ed is simply given by X ˆ S:

…6†

This condition is equivalent to vanishing distance between the saddle S and the point X on the axis (D). Mathematically speaking, conditions (4) and (6) are equivalent. Note, however, that condition (6) is accessible via approximations of periodic orbits. The Melnikov condition (4), on the other hand, circumvents periodic orbits entirely, aiming at the separatrices directly. 3. Elliptic averaging method Consider Eq. (2), where dots denote derivatives with respect to time, f …x† ˆ x2 ; g…l; x; x_ † is an arbitrary nonlinear function of its arguments and l is referred to as a control parameter. A survey of elliptic function properties is given in Appendix A. Eq. (2) has two equilibria, …0; 0† and …ÿc1 =c2 ; 0† whose stability depends on c1 and c2 . To apply the elliptic averaging method [15] to Eq. (2), we ®rst solve the unperturbed system (3) with the Ansatz 2

x ˆ A1 ‡ A2 …sn† ;

u ˆ xt ‡ b;

…7†

in which sn ˆ sn…u; m† and A1 ; A2 ; x; b and m are constants. Substituting Eq. (7) and its derivatives into Eq. (3) yields 2

4

S1 ‡ S2 …sn† ‡ S3 …sn† ˆ 0;

…8†

where S1 ˆ 2x2 A2 ‡ c1 A1 ‡ c2 A21 ; S2 ˆ A2 ‰ÿ4x2 …1 ‡ m† ‡ c1 ‡ 2c2 A1 Š;

…9†

2

S3 ˆ A2 …6mx ‡ c2 A2 †: Requiring S1 ; S2 and S3 to vanish gives three nonlinear algebraic equations relating the six unknown parameters A1 ; A2 ; x; b; m and l …c1 and c2 are assumed to be known). Solving for A1 ; A2 and x in terms of m; c1 and c2 provides

2254

M. Belhaq, F. Lakrad / Chaos, Solitons and Fractals 11 (2000) 2251±2258



1=4 c21 xˆ ; 16k p p ÿc1 k ‡ c21 …1 ‡ m† p ; A1 ˆ 2c2 k r ÿ3m c21 ; A2 ˆ k 2c2

…10† …11† …12†

where k ˆ 1 ÿ m ‡ m2 . The other three unknown relations are given by the modulation equations of m and b and the collision criterion (6). On the other hand, note that the apparently more general system (2) can be reduced by ane transformations on x and t, to the following equation: x ‡ x ‡ x2 ˆ eg…l; x; x_ †;

…13†

which will be considered in what follows. Since the unperturbed system is conservative, it admits the energy integral. x_ 2 x2 x3 ‡ ‡ ˆ h: 2 2 3

…14†

After substituting the expressions of x and its derivative in Eq. (14), we ®nd the following relation between energy h and square-modulus m:   1 …m ÿ 2†…m ‡ 1†…2m ÿ 1† 1ÿ : …15† hˆ 12 2k3=2 Following [15] and using the variation of parameters to express the e€ect of the order e terms on the slow evolution of the square-modulus m, now considered a function of t, we ®nd the equation m_ ˆ

eg…l; x; x_ †xb ; xmt xb ÿ xbt xm

…16†

in which subscripts represent partial di€erentiation. Macsyma is used [15] to do substitutions and simpli®cations of various identities to ®nd m_ ˆ e

8 k9=4 x_ g…l; x; x_ †: 9mx …1 ÿ m†

Application of the method of averaging to (17) yields Z 4K 8 k9=4 1 x_ …u†g…l; x…u†; x_ …u†† du: m_ ˆ e 9mx …1 ÿ m† 4K 0

…17†

…18†

An equilibrium point of the averaged equation (18) corresponds to a limit cycle in the original equation (13). Thus if m ˆ m0 is a root of the equation Z 4K x_ …u†:g…l; x…u†; x_ …u†† du ˆ 0; …19† 0

the averaged equation predicts that for small e, a limit cycle coincides with the energy curve, associated with a value of h ˆ h0 which corresponds to m0 . The coordinates of the saddle in Eq. (2) depend on the choices of c1 and c2 . Two cases are possible Case 1: c1 < 0 and hence …as ; bs † ˆ …0; 0†: In this case, the homoclinicity criterion (6) can be written as x…u† ˆ 0; x_ …u† ˆ 0:

…20† …21†

M. Belhaq, F. Lakrad / Chaos, Solitons and Fractals 11 (2000) 2251±2258

2255

To leading order, these equations are given explicitly by the system A1 ‡ A2 …sn†2 ˆ 0;

…22†

2A2 …sn†…cn†…dn† ˆ 0:

…23†

To ensure A1 6ˆ 0 and A2 6ˆ 0, the condition sn 6ˆ 0 is required. In order to satisfy Eq. (23) we must impose either the condition cn ˆ 0 or dn ˆ 0. The ®rst condition leads to A1 ˆ ÿA2 ;

…24†

and requires m ˆ 1. The second condition yields A1 ˆ ÿ

A2 ; m

…25†

which requires also m ˆ 1. Case 2: c1 > 0 and hence …as ; bs † ˆ …ÿc1 =c2 ; 0†: In this case the homoclinicity criterion (6) reads c1 ; c2 2A2 …sn†…cn†…dn† ˆ 0: 2

A1 ‡ A2 …sn† ˆ ÿ

…26† …27†

In order to satisfy A1 6ˆ 0 and A2 6ˆ 0 it is required that cn ˆ 0 or dn ˆ 0. Both cases lead to m ˆ 1 as before. Hence, the homoclinicity condition (6) is equivalent to m ˆ 1. For an arbitrary polynomial g…l; x; x_ †, the limit cycle integral condition (19) may be eciently evaluated by using computer algebra. However, the integration from 0 to 4K in (19) can be shifted to the interval t 2 ‰ÿ2K; 2KŠ because the integrand is of period 4K in u. This indeed produces the integral over u 2 …ÿ1; ‡1† in the limit K % 1 which corresponds to the collision between the saddle and the limit cycle given by m % 1. Consequently, Eq. (19) obtained via averaging coincides with the Melnikov function (see [25]). 4. Applications As a ®rst example, we consider the quadratic Arnold±Bogdanov±Takens equation [29] u_ ˆ v; v_ ˆ l1 ‡ …l2 ‡ u†v; p with l1 < 0. Setting x ˆ u ‡ ÿl1 Eq. (28) becomes p p x ‡ 2 ÿl1 x ÿ x2 ˆ e…l2 ÿ ÿl1 ‡ x†_x:

…28†

…29†

The small parameter e is introduced to have the form of Eq. (2). The ®xed points of (29) are a focus …0; 0† p p p and a saddle …2 ÿl1 ; 0†. Here c1 ˆ 2 ÿl1 ; c2 ˆ ÿ1 and g…l; x; x_ † ˆ …l2 ÿ ÿl1 ‡ x†_x: The solution is sought in the following form: x ˆ A1 ‡ A2 sn2 …xt ‡ b; m†: Then, Eqs. (10)±(12) lead to   …1 ‡ m† p ÿl1 ; A1 ˆ 1 ÿ p k p 3m ÿl1 p ; A2 ˆ k  ÿl 1=4 1 : xˆ 4k

…30†

…31† …32† …33†

2256

M. Belhaq, F. Lakrad / Chaos, Solitons and Fractals 11 (2000) 2251±2258

To determine l2 ; m and b we use the conditions of stationarity of m and b and the condition of homoclinicty given by m % 1. Hence, Eq. (19) leads to the relation p C2 ‰4…1 ÿ m†A22 …l2 ÿ ÿl1 † ‡ 4…1 ÿ m†…A1 ‡ A2 †A22 Š p ‡ C4 ‰4A22 …l2 ÿ ÿl1 †…2m ÿ 1† ‡ 4…A1 ‡ A2 †A22 …2m ÿ 1† ÿ 4A32 …1 ÿ m†Š p …34† ‡ C6 ‰ÿ4A22 m…l2 ÿ ÿl1 † ÿ 4m…A1 ‡ A2 †A22 ÿ 4A32 …ÿ1 ‡ 2m†Š ‡ C8 4mA32 ˆ 0; which gives when m % 1 A1 ‡ …l2 ÿ

3 p ÿl1 † ‡ A2 ˆ 0: 7

…35†

Substituting the expressions of A1 and A2 , given in Eqs. (31) and (32), into (35) provides the homoclinic bifurcation curve  2 7 l22 : …36† l1 ˆ ÿ 5 This prediction coincides with the result given by the Melnikov technique [29]. In the second example we consider the generalized van der Pol oscillator x ‡ x ÿ x2 ˆ e…l ÿ x2 †_x:

…37†

The ®xed points of Eq. (37) are a focus …0; 0† and a saddle …1; 0†. Here c1 ˆ 1; c2 ˆ ÿ1 and g…l; x; x_ † ˆ …l ÿ x2 †_x. Therefore, condition (19) becomes 2

2

C2 ‰4A22 …1 ÿ m†…l ÿ …A1 ‡ A2 † †Š ‡ C4 ‰4A22 …ÿ…2m ÿ 1†…A1 ‡ A2 † ‡ 2A2 …A1 ‡ A2 †…1 ÿ m† ‡ l…2m ÿ 1††Š ‡ C6 ‰ÿ4A22 …lm ‡ …1 ÿ m†A22 ÿ 2A2 …2m ÿ 1†…A1 ‡ A2 † ÿ m…A1 ‡ A2 †2 †Š ‡ C8 ‰ÿ4A22 ……2m ÿ 1†A22 ‡ 2mA2 …A1 ‡ A2 ††Š ‡ C10 4mA42 ˆ 0:

…38†

In the limit m % 1 of Eqs. (38), (10)±(12) we obtain the following critical value of homoclinic bifurcation: 1 lˆ ; 7

…39†

which coincides also with Melnikov literature [30]. 5. Conclusion In a recent paper [22], a new analytical technique to derive a criterion for predicting homoclinic bifurcations in autonomous dynamical systems was presented. This criterion is mainly based on the collision between the bifurcating periodic solution and the saddle equilibirium. In this approach, however, the periodic orbit was approximated using trigonometric perturbation methods. Due to the local character, near the ®xed point of the trigonometric approximations to the periodic orbit, the obtained results furnished reasonable predictions of homocline bifurcations. In this work, we have combined the collision criterion with the Jacobian elliptic functions to establish a homoclinic bifurcation criterion taking advantage of the global character of elliptic functions. As an asymptotic expansion we have adopted the averaging elliptic method. We have shown that the results obtained by this perturbation technique agree with those predicted by the Melnikov method to leading order. A similar approach combined the collision criterion and the Jacobian elliptic Lindstedt-Poincare perturbation method was conducted to predict homoclinicity in autonomous systems [25]. The results of this technique coincide also with the Melnikov function. Despite the slight di€erences in the two analytical perturbation methods, the elliptic Lindstedt±Poincare and the elliptic averaging, the obtained results are the same.

M. Belhaq, F. Lakrad / Chaos, Solitons and Fractals 11 (2000) 2251±2258

2257

A natural question raised here is how one can take advantage of these alternative techniques to derive a second order approximation to homoclinic bifurcation. Extensions of the method to predict homoclinic bifurcations in three-dimensional systems is also in order. A ®rst tentative step in this direction will be published elsewhere. Appendix A For the convenience of our readers, we collect some facts on Jacobian elliptic functions [31]. Jacobian elliptic functions satisfy algebraic relations which are analogous to those for trigonometric functions. The fundamental three elliptic functions are cn…u; m†; sn…u; m† and dn…u; m†. Each of the elliptic functions depends on the square of the modulus m as well as the argument u. Note that the elliptic functions sn and cn may be thought of as generalizations of sin and cos where their period depends on the modulus m. The elliptic functions satisfy the following identities, which are analogous to sin2 ‡ cos2 ˆ 1: sn2 ‡ cn2 ˆ 1; msn2 ‡ dn2 ˆ 1; 2

…A:1† 2

1 ÿ m ‡ mcn ˆ dn : Only two of these relations are algebraically independent. In Table 1, additional properties of Jacobi elliptic functions are summarized. Here K…m† is the complete elliptic integral of the ®rst kind K…0† ˆ p=2;

K…1† ˆ ‡1:

We de®ne Z C2P ˆ

4K…m†

0

cn2P du:

…A:2†

Then C0 ˆ 4K…m†;

C2 ˆ

4 ‰E…m† ÿ …1 ÿ m†K…m†Š: m

…A:3†

Here E…m† is the complete elliptic integral of the second kind. In general, we have the recursion relation C2P ‡2 ˆ

2P 2m ÿ 1 2P ÿ 1 1 ÿ m C2P ‡ C2P ÿ2 : 2P ‡ 1 m 2P ‡ 1 m

…A:4†

When m increases from 0 to 1, E…m† decreases from p=2 to 1. In the limit m % 1, the quantities K…m† and E…m† have the following behavior, to leading order: K…m† ˆ

1 16 log ‡


; 2 1ÿm

E…m† ˆ 1 ‡




…A:5†

Table 1 Additional properties of Jacobi elliptic functions Property

sn…
; m†

sin…
†

cn…
; m†

cos…
†

dn…
; m†

Max. value Min. value Period Parity df =du fmˆ0

1 ÿ1 4K…m† Odd cn
dn sin

1 ÿ1 2p Odd cos sin

1 ÿ1 4K…m† Even ÿsn
dn cos

1 ÿ1 2p Even ÿ sin cos

1  p 1ÿm 2K…m† Even ÿm
sn
cn 1

2258

M. Belhaq, F. Lakrad / Chaos, Solitons and Fractals 11 (2000) 2251±2258

In particular, 4K…m†, the minimal period in integral (19), approaches in®nity in the homoclinic limit m % 1. With recursion (A.3) and (A.4), the asymptotics (A.5) imply lim C0 ˆ 1;

m!1

lim C2 ˆ 4;

m!1

8 lim C4 ˆ ; 3

m!1

lim C6 ˆ

m!1

32 ; 15

lim C8 ˆ

m!1

192 ; 105

lim C10 ˆ

m!1

512 ; 315

…A:6†

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31]

Nayfeh AH, Mook DT. Nonlinear oscillations. New York: Wiley, 1979. Jordan DW, Smith P. Nonlinear oordinary di€erential equations. Oxford: Oxford University Press, 1987. Nayfeh AH. Perturbation methods. New York: Wiley, 1973. Nayfeh AH. Introduction to perturbation techniques. New York: Wiley, 1981. Krylov N, Bogolioubov N. Introduction to nonlinear mechanics. Princeton, NJ: Princeton University Press, 1943. Bogolioubov N, Mitropolsky I. Asymptotic methods in the theory of nonlinear oscillations. New York: Gordon and Breach, 1961. Kevorkian J, Cole JD. Perturbation methods in apllied mathematics. New York: Springer, 1981. Barkham PGD, Soudack AC. An extension to the method of Krylov and Bogolioubov. Internat J Control 1969;10:377±92. Barkham PGD, Soudack AC. Approximate solutions of nonlinear non-autonomous second-order di€erential equations. Internat J Control 1970;11:101±14. Soudack AC, Barkham PGD. Further results on approximate solutions of nonlinear, non-autonomous second-order di€erential equations. Internat J Control 1970;12:763±7. Soudack AC, Barkham PGD. On the transient solution of the unforced Dung equation with large damping. Internat J Control 1971;13:767±9. Yuste SB, Bejarano JD. Extension and improvement to the Krylov±Bogolioubov methods using elliptic functions. Internat J Control 1989;49:1127±41. Bejarano JD, Sanchez AM, Rodriguez CM. Osciladores alineales excitados no linealmente. Annales de Fisica A 1982;78:159±64. Yuste SB, Bejarano JD. Amplitude decay of damped nonlinear oscillators studied with Jacobian elliptic functions. J Sound Vibration 1987;114:33±44. Rand RH. Using computer algebra to handle elliptic functions in the method of averaging. Symbolic computations and their impact on dynamics. American Soc Mech Engrg PVP 1990;205. Yuste SB, Bejarano JD. Construction of approximate analytical solutions to a new class of nonlinear oscillator equations. J Sound Vibration 1986;110:347±50. Garcia-Margallo J, Bejarano JD. Generalized Fourier series and limit cycles of generalized van der Pol oscillators. Internat J Control 1988;49:1127±41. Coppola VT, Rand RH. Averaging using elliptic functions: approximation of limit cycles. Acta Mech 1990;81:125±42. Coppola VT, Rand RH. Macsyma program to implement averaging using elliptic functions. Computer aided proofs in analysis. Germany: Springer, 1991;71±89. Chen SH, Cheung YK. An Elliptic Lindstedt±Poincare method for certain strongly nonlinear oscillators. Nonlinear Dynamics 1997;12:199±213. Chen SH, Cheung YK. An elliptic perturbation method for certain strongly nonlinear oscillators. J Sound Vibration 1996;192(2):453±64. Belhaq M. New Analytical technique for predicting homoclinic bifurcations in autonomous dynamical systems. Mech Res Comm 1998;23(4):381±6. Belhaq M, Fahsi A, Lakrad F. Predicting homoclinic bifurcations in planar autonomous systems. Nonlinear Dynamics 1999;18(4):303±10. Belhaq M, Fahsi A. Homoclinic bifurcations in self-excited oscillators. Mech Res Comm 1996;23(4):381±6. Belhaq M, Fiedler B, Lakrad F. Homoclinic connections in strongly self-excited nonlinear oscillators: the Melnikov function on the elliptic Lindstedt±Poincare method, Nonlinear dynamics, to appear. Belhaq M, Houssni M, Freire E, Rodriguez-Luis AJ. Asymptotics of homoclinic bifurcation in a three-dimensional system. Nonlinear Dynamics, to appear. Chow N, Hale JK. Methods of bifurcation theory. New York: Springer, 1982. Kuznetsov YA. Elements of applied bifurcation theory. New York: Springer, 1995. Guckenheimer J, Holmes PJ. Nonlinear oscillations, dynamical systems and bifurcation of vector ®elds. New York: Springer, 1983. Merkin JH, Needham DJ. On in®nite period bifurcations with an application to roll waves. Acta Mech 1986;60:1±16. Byrd P, Friedman M. Handbook of elliptic integrals for engineers and scientists. Berlin: Springer, 1971.