The approach to constrained equations

Physics Letters A 183 ( 1993 ) 319-326 North-Holland

P H Y S I C S IIZTTI:RS A

The approach to constrained equations Darryl Veitch Telecom Research Laboratories, 770 Blackburn Road, Clayton, Victoria 3168, Australia Received 11 May 1993; accepted for publication 12 October 1993 Communicated by A.R. Bishop

We consider the singular perturbation of a constrained equation in three dimensions with a two-dimensional "slow" manifold, to a three-dimensional differential system with two time scales. Via a general methodology involving a reduction to one-dimensional maps, we obtain numerically a growth law which describes the transition quantitatively.

1. Introduction Apart from its intrinsic interest, the theory of constrained equations has received considerable attention from its potential to explain the relaxation type oscillations [1], present in many physical systems which possess two widely disparate time scales. Such systems can be described by equations of the form ki = f ( x l ..... x,,), Yci= -1 hi(x1 ..... x , ) , E

i=l,...,k,

i = k + l , ..., n ,

(1)

where E is the ratio of the slow to the fast time scale. In contrast, constrained equations are not differential equations but dynamical systems of the following form, . ~ i = f i ( X l ..... Xn),

i=1 ..... k,

0=h~(x~,...,xn),

i = k + 1 ..... n .

(2)

Through the implicit function theorem, there are local existence theorems for such equations, but global solutions require special attention. In the traditional picture, which Takens formalised in his works [ 2,3 ], global solutions consist of differentiable flows on the stable branches of the submanifold M defined by the constraints, with non-differentiable, instantaneous jumps between different branches. Except possibly initially, these jumps occur only at singularities where the stable and unstable branches meet. The "slow" manifold of eqs. (1), where the fast variables ,are

slaved to the slow, and which is approximated by M for small E, is replaced by M in the limit. Global solutions defined in this way were taken to meaningfully model those of physical systems with two different time scales in the singular limit e-,0. The assumption that this is valid, sometimes referred to as the "perfect delay convention" [4 ], is the classical one made by RSssler in his models of abstract chemical reaction kinetics [ 5 ], and electrical circuits [ 6 ], and by Tyson [ 7 ] and Rinzel and Troy [ 8 ] (to name a very few), in their studies of the BelousovZhabotinskii reaction. It was soon realised however, that for arbitrarily small E, systems with two time scales could display complicated solutions, which did not approach any member of the class of intuitive solutions of constrained equations just described. In refs. [9,10] non-standard analysis ,1 is employed to demonstrate the existence of and describe these non-intuitive "canard" solutions. Canards are solutions which travel over both the stable and unstable branches of the slow manifold. Their existence necessitates a much broader definition for solutions of constrained equations, one without uniqueness of trajectories and where the unstable branches play a full role. In refs. [ 13,14 ] Benoit provides a comprehensive analysis of the singular perturbations of typical three-dimensional constrained equations (with ~1 Wewillnot actuallyemploynon-standardanalysisas suchhere. The reader is referred to refs. [ 11,12,10] for information on this fascinatingsubject.

0375-9601/93/$ 06.00 © 1993 Elsevier Science Publishers B.V. All rights reserved.

319

Volume 183, number 4

PHYSICS LETTERS A

two-dimensional slow m a n i f o l d s ) . He shows that in the great m a j o r i t y o f cases the relationship between the solutions o f the perturbed systems (which are o f the form o f eqs. ( 1 ) ), and the b r o a d e r class o f constrained equation solutions designed to describe them, is one to one. It has been noted however [9 ], that canards are very difficult to observe, and are very likely to be sensitive to noise [4]. O f course in addition, and as we shall see, they will not always be relevant to the d y n a m i c s of interest. Thus for m a n y practical purposes one m a y be able to think in terms o f the intuitive solutions only. In this paper we deal quantitatively with the relationship between constrained systems and their singular perturbations. We a d o p t a numerical approach, studying the asymptotic dependence o f the fast-slow flow on e in a typical family o f threed i m e n s i o n a l constrained equations with a twod i m e n s i o n a l slow manifold, M. In order to do this we need a way to measure quantitatively the variation o f solutions o f the p e r t u r b e d equations with E. This is achieved by characterising the family o f flows ( p a r a m e t e r i s e d by E) by a corresponding family of one-dimensional m a p s #2 which accurately describe much o f their dynamics. Certain stable periodic orbits (trajectories) o f these m a p s are the specific objects which we finally quantitatively compare. O u r m a i n result is a growth law describing how solutions m o v e away from their constrained equations' limit after perturbation. The form o f the law and its p a r a m e t e r (the growth rate) are characteristic o f the particular constrained equation which underlie the physical system. The methodology presented is generally applicable to other constrained systems and suggests the possibility o f a classification o f constrained systems by growth rate. In section 2 we construct an explicit family o f constrained equations which possess interesting return ~2 These one-dimensional maps are return maps, defined by the action of trajectories on M as they cross a suitable linear section. This method of effectively reducing the dimensionality of a system is standard practice in dynamical systems theory, and has been applied in the context of constrained equations many times (for example see refs. [ 15,3,16,17,8 ] ). In the case of the perturbed equations, a two-dimensional return map can be defined transversal to M which is "quasi"-one-dimensional, because of the strong compression of the flow onto the slow manifold. We form a one-dimensional map from it simply by discarding the fast coordinate. 320

13 December 1993

maps, in section 3 we integrate them numerically and locate interesting p a r a m e t e r regimes, a n d in section 4 we singularly perturb the constrained system and calculate the growth law.

2. The constrained equations To obtain interesting dynamics, we need to design a system with a b r o a d m e c h a n i s m for recurrence. This is not difficult to achieve if we have the freed o m to postulate constraining surfaces and the flows on t h e m at will, a n d m a n y interesting return maps describing such systems have been generated in this way (for example, see refs. [ 3 , 1 7 ] ) . It is another matter however to construct explicit equations to m o d e l an envisaged mechanism. Consider the following constrained equations. They i m p l e m e n t a surface/flow combination, illustrated in fig. 1, which possesses recurrent dynamics described by a b i m o d a l return m a p (one with two turning points).

i = a ( z ) [X-Xo(Z) ] + w[y-yo(Z) l = f ( x , y, z) ,

(3a)

Fig. 1. The basic flow/fold combination embodied in the constrained equations ( 3a )- (3c). The section Zt is mapped into itself after jumping from Mt+ to M~ and returning.

Volume 183, number 4

PHYSICS LETTERS A

13 December 1993 F6:

S'= - w [ x - x o ( z ) ] +a(z) [y-yo(Z) ]

[Ft .

(3b)

0=h-(z) -x=h(x, z),

,

. •.

..

".

(3c)

i,,'

where h-(z) = - a z ( z + d ) ( z - d )

,

Xo(Z) = j ( z ) (Xb --Xt) +Xt, yo(Z)=j(z)yb,

.....

•..

,

=g(x,y, z ) ,

.•

a(z)=j(z)(ab--Oq)+Ott,

Ct

' C~

(a)

j(Z) = arctan( - zsz) /lt + ½ and where w, a, d, s=, ab, at~R + and xb, xt, yb~R. Note that: (i) f g and h are C 1 and h = 0 defines a smooth two-dimensional manifold, (ii) the implicit setting of y, to zero incurs no loss of generality, (iii) j ( z ) maps ( - o o , oo) to (1, 0) and (iv) there are three singular points: At at (0, 0, d), Ab at (0, 0, - d ) , and the origin. Let M + = M ~ - w M +, where M + ={{x, y, z}eM: z>dx/~ } denotes the top stable branch of M, and M~- the bottom stable branch ( z < - d v ~ ) . The spiral trajectory shown on M + in fig. 1 touches the foldline F~ at I-It, a pseudo-singular point of saddle type. (We will ignore the implications of this for the present, see ref. [ 13 ].) Hence it separates the trajectories which have infinite histories on M + from those which have not, and thus provides an optimal choice of linear section, Zt. If this "limiting" spiral trajectory, and the analogous one for M~-, are positioned as shown, then an invariant region is created on M + within which all trajectories except At and Ab cycle endlessly between M + and M~-, and Zt is repeatedly mapped back into itself• This can be seen in projection in fig. 2, where the limiting spirals and the return map are plotted from numerical solutions of eqs. (1). For a suitably large s=, j(z) is approximately zero on Mt+ and unity on M ~ , and so (Xo, Yo, a ) ~ (xt, Yt, at) on M + , and (xb, Yb, orb) on M~-. In this way the horizontal components (f(x, y, z), g(x, y, z) ) are distinct and conceptually simple on each branch of M +, yet are smoothly joined, so that the singularly perturbed form of the system, formed by replacing eq. (3c) by k=h(x, z)/e, is a well defined set o f ordinary differential equations in three dimensions. Arbitrary flows could be joined in this

..... I.... Cb

Ebi

.-'•

(b)

Fig. 2. Numerical calculations of eqs. (3a)-(3c) with (X't, at, X{,, ab) = ( 1.0, 0.7, - 1.0, 0.7). (a) The projections of the foldlines, the sections and the limiting spirals, and the trajectories which define the turning points Ct and Cb are shown. Dashed lines indicate the lower surface• (h) The return map R over the section Zt (250 points)• It possesses a stable periodic orbit and three turning points: Ct and two preimages of Cb.

321

Volume 183, number 4

PHYSICS LETTERS A

manner, but we only vary coefficients o f a topologically unchanged flow, that o f a linear planar (unstable) focus. Thus we choose a second spiral on M~- to generate the desired reinjection to M~+ , rather than, for instance, a unidirectional flow as in refs. [ 18,17 ]. We do this for two reasons. First, this configuration generates b i m o d a l dynamics o n Et, because the maps Y ~ t ~ b and Y.b~Et each have one fold, due to the fact that the projection o f F t ( F b ) is tangent to a unique trajectory on M~- ( M + ) . See the points C~ (C{,) in fig. 2. Bimodal m a p s ~3 are o f considerable interest because they are non-trivial and yet well u n d e r s t o o d d y n a m i c a l systems. They can possess two i n d e p e n d e n t attractors, a feature o f great i m p o r t a n c e in the observation o f physical systems. They are essential to the understanding o f non-invertible circle maps, which have been used extensively to model the transitions to chaos in two-frequency systems [20]. The second reason for only varying coefficients is that, since it is possible to generate interesting dynamics with only a modest change in the horizontal flow, it is more natural to do so than to choose to "forcibly" j o i n two very different ones. U n d e r more general conditions a more controlled evolution o f the horizontal c o m p o n e n t o f the field with z is not just desirable, but m a y be necessary for recurrence. F o r instance in the case when the return to M t is achieved by means o f a flow a r o u n d a cusp rather than by another j u m p . The theory o f piecewise continuous m a p s on the interval and b i m o d a l maps in particular, is quite well developed. It is well known, see refs. [21,22], that the dynamics up to topological conjugacy is determ i n e d by the " k n e a d i n g invariant", that is the forward trajectories o f the turning points. Ideally then, one would wish to establish the relationship between the p a r a m e t e r s in eqs. ( 3 a ) - ( 3 c ) a n d these two internal parameters. This is no simple task however, so we turn to numerics to generate return m a p s and exp e r i m e n t to find interesting p a r a m e t e r regions. Note that the return m a p R: gt--,Zb in fig. 2b, has three turning points: Ct is the first preimage o n 'Et o f Ct, and P~ (Pr), are the left (right) first preimages on E t o f C b ~ E b. This is a general feature, since R is de~3 This work formed part of a doctoral thesis [ 19 ] which dealt mainly with bimodal maps in various contexts. 322

I

13 December 1993

[ C,

[

E,

I Cb

I

r.b

Fig. 3. The graphs of the maps gt~Y-~ and Y-q,~Ztwhich define the underlying bimodal map. The "second iterate" mappings Z ~ Z i are also shown. We denote by R the case with i= t and deal with it exclusively. fined over an invariant interval o f the second iterate o f the underlying b i m o d a l m a p shown in fig. 3.

3. N u m e r i c a l integration of the constrained equations

The constrained equations ( 3 a ) - (3c) were integrated as a set o f second order differential equations on M + with instantaneous j u m p s between the stable surfaces inserted at the folds. The integration was performed in (z', x) space rather than in (x, y ) space, where z'= ( z - F ~ ) 2 is a quadratic z variable centred on F~, thus avoiding the necessity o f solving a cubic equation at each t i m e step for z. All numerical work was carried out on an IBM 3084 m a i n f r a m e computer a n d was written in Fortran. There are m a n y parameters which could be used to control R. We fix w = 3, a = 5, d = 2 , and set Sz equal to 100 times the vertical distance o f the origin to the folds. Further, we make Yb a function o f the other parameters by setting it such that lib coincides (in projection) with the image o f lit on M + as shown in

Volume 183, number 4

PHYSICS LETTERS A

13 December 1993

fig. 2a. This choice was motivated by numerical experiments and helps ensure that Et is invariant, as the limiting spirals do not fit as nicely as in fig. 2 for all parameter values. Thus we write R (c~t, ab, Xt, Xb) : Et"+ E b .

We search (at, ab, xt, Xb) space for interesting parameter regions, that is where R (c~t, Orb, Xt, Xb) possesses a stable periodic orbit o f fairly low period (but greater than 2). Not only are such orbits interesting from a physical point of view, relatively easy to find numerically, and unlikely to be artifacts o f the numerics, but they also give us an easy way of determining the orbits of the turning point(s) involved. This is because turning points are always (with unimportant exceptions) attracted to stable periodic orbits, if they exist. O f particular interest are regions in parameter space where two stable periodic orbits exist simultaneously, each attracting the orbit of a different turning point. Especially if they have different periods, this behaviour exposes very unambiguously that the underlying behaviour is due to a map which is at least bimodal. Experimentally if such a region was found in a relaxation oscillator it would be very suggestive of underlying bimodal dynamics. It is generally the case that there is an invariant attracting subinterval in the interior of Et. Canard solutions, if they ever enter this subinterval, can never escape and will therefore never return to the neighbourhood of Fit or rib, the pseudo-singular points. On the other hand if they never enter, they must be unstable. Thus for the parameter regimes we consider, one can safely ignore canards from the point of view of observable dynamics. In this system the regions where stable periodic orbits exist are very narrow, as the slopes of R are large away from the turning points. Hence regions where two coexist are very small indeed (like a small open neighbourhood of a co-dimension two set). We locate points in these regions using an algorithm which searches in (At, C~b) space, given only the orders of the periodic orbits desired and values of (at, Xb) which enable the search to succeed (determined by experiment). Two examples are given in figs. 4 and 5. The first few iterates of each turning point are shown, highlighted with "graphical analysis" lines. In addition, for each periodic orbit the first 100 iterates of the associated turning point are shown as dots below the

Fig. 4. The return map R(X't, at, X{,, Ctb, Y{,)=R(8.96523, 2.0, 9.0, 1.56320, 95.8700) at a point in parameter space where stable periodic orbits of orders 3 (Ct) and 5 (Pr) coexist.

"/

:.

R"

P Fig. 5. The return map R (X't, at, X{,, ab, Y{,)= R ( - 7.7979, 1.2, 8.0, 1.3284, 128.73) at a point in parameter space where each of Ct and P~ control dynamically distinct unimodal maps, which possess stable periodic orbits of period 3 and 4, respectively. 323

Volume 183, number 4

PHYSICS LETTERS A

return map. The periodic orbits are very close to the orbits of the turning points, because the former have very small immediate basins of attraction. In fig. 4 the dynamics is truly bimodal, as the periodic orbits are linked, whereas in fig. 5 it reduces to that of twounimodal maps. The algorithm searches for a periodic orbit of period q which attracts a turning point C, by looking for an approximation to a periodic orbit which actually passes through C. Such orbits are known as "super-stable". When a good enough approximation is found the map will possess a stable periodic orbit of the desired order which passes near to C, and which attracts the orbit of C. If we wish to actually locate the super-stable point we need only increase the accuracy demanded of the approximation. To carry out the search, C - R q ( C ) is passed to a root solving algorithm acting in a one-dimensional parameter space. The parameter C~b was used to search for orbits around C~, as experiment had shown a reasonably linear connection between them in this regime. For the search based on P~ or Pr the parameter xt was used. This combination was found to be orthogonal enough so that it was possible to employ the two one-dimensional searches alternatively until convergence for both simultaneously was achieved (a direct two-dimensional root solving algorithm was tried but had poor convergence properties). Locating these points both in parameter space ( ~ sixth digit) and in phase space ( ~ fifth digit) was a delicate affair. It was therefore essential to develop an efficient method of locating the desired orbits, so that it would be feasible to pursue them into ~ space.

13 December 1993

unfortunately completely mistaken, as the slaving involves a delicate balance where the distance of the trajectory to M + is critical. It is therefore necessary at all times to integrate the fast component accurately, which with ordinary schemes results in an extremely slow integration over M +. Fortunately the state of the art numerical scheme for dealing with such systems is encoded in the NAG subroutine library stiff-equation routines. These indispensable tools make it possible to integrate the system to high accuracy ( ~ 8 digits) and for very high stiffness ratios, in fact as high as E = 10 - 9 . To give an idea of the variation with ~ of the return map at fixed parameter values we present fig. 6. It shows graphs of R(8.9, 2.0, 9.0, 1.495) for the constrained equation, and for the differential equation with e=0.001 and 0.01. For the constrained equation, R possesses a stable periodic orbit of order three passing near Ct which does not persist even for the smaller perturbation: e = 0.001. To investigate this in detail we take the opposite approach to fig. 6 and vary parameters to maintain the behaviour of R, rather than leave the parameters unchanged and try to measure how the dynamical behaviour varies. To do this

-.. ".:'.

,/" j/

..;.. "i:.

,

""

4. The singular perturbation Computationally the integration of the constrained and perturbed equations are very different, despite the convergence in their behaviour as E--.0. By definition the differential system is "stiff", that is it possesses widely varying time scales. Such systems are notoriously difficult to integrate numerically with any accuracy. Sophisticated specialised techniques are required, for example the Gear algorithm [23]. One might think naively that careful integration would only be required near the folds, as the fast motion is slowed to the slow on M ÷. This is 324

Fig. 6. Three return maps for different ~, with (X;, o~t, X~,, O~b)= (8.9, 2.0, 8.5, 1.495). The denser plot is the constrained equation, then e=0.001 (close to the dense plot), and finally e=O.Ol.

Volume 183, number 4

PHYSICS LETTERS A

qualitatively we must focus on a sharper feature than just a stable periodic orbit, as they exist for open parameter regions. By insisting on super-stability of both orbits we specify the maps very precisely, in fact generically to just a point in a two-parameter search. We only need however to follow a single super-stable orbit in a oneparameter search to be equally specific. Three searches in ab (E) were made for orbits including Ct

g.SL~l*

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

i,

,

, + , ,

13 December 1993

of periods 3, 4, and 5. Figure 7ais an ( a b ( e ) , e) plot for the period 3 case (the others look similar), the log-log plots then follow in order in figs. 7c, 7d. For these latter plots the data has been translated to make otbo (the ab value of the constrained equation) the origin. The range o f e over which the points were calculated was e = 10-9-10-3• In each case the log-log graphs were very close to straight lines and the slopes of each, calculated from the five points with smallest

......

I.E-4~

J

l

........

.......

J

.......

.I

,

.......

l

.......

d

/ II I I 1 5

(a)

I MII

1 .E.,416

II.IIIS

1 t..4m

I.lUl

I I.II

1.21

1.22

(b)

I ,IE-ll

1-24

1.16

I•~II

1.1

........ I

..... I

i.~l

~..IEll

....... I

I•~"14

~..IE45

~b

........

1.01[--413

I

.......

I

........

I

........

-1

.....

,.I

1 .IE--II4

I .IE-0S

....... I l,IEll

......

I. I K 4 3

J

........

I

.......

....... I ........ I .IE-4111 I •IE-~411

(O~b - - ~ b 0 ) J

........

I

........

I

I .I~-14 •

•

I IE,4ll

I 0E-II6

I E-It

I 0£.-ll7

I I~,-17

(c)

(d) 1 .IE-II

-

I• E l i

1.0E-09

I.IE-ll I. I E l l

-

l.IK-l~

........ I I,I£-04

.......

I

I .IE-IIJ

........

I

l ,lI:'-4~I

........

I

I .E-Ill

........

I .E-ell

I •IE-II l.Ell

--(Or b --Olb0 )

......

l.E--If

"]

........

I

I. I E 4 I

........

I.

1 E l l

........

I

.......

l, I ~ 4 1

I JK.t411

- - (O~b - - ~ b 0 )

Fig. 7. The results of following super-stable orbits in ~ space. The searching variable is ~b in each case and C~ is the turning point which is periodic. (a) Plot of the actual values of ~versus ab for a period 3 example; (b) the corresponding log-log plot; (c), (d) log-log plots

for examples with orders 4 and 5, respectively. The slope in each case is 1•50• 325

Volume 183, number 4

PHYSICS LETTERSA

smallest E to obtain the limiting behaviour, was 1.50 with uncertainty in the third decimal place. This strongly suggests a dependence of the form -

( a b - a b o ) ~ e 2/3 •

(4)

Our claim is that this relation, which is characteristic of the change in R with e << 1 (not just of certain trajectories) and hence of the change in solutions of the perturbed system with ¢, is a function of the underlying constrained system and could be used as an identifying signature for it. 5. Conclusion

We have shown that a particular system with two time scales is characterized by a growth law ~2/3 where e measures the distance from the limiting constrained system. Assume for the m o m e n t that different classes of constrained systems have different growth law signatures or growth "rates". Now if e is an accessible parameter in systems displaying relaxation oscillations (which seems a reasonable ass u m p t i o n ) , such a quantity could be measured experimentally provided (for instance) that suitable stable periodic solutions could be located. If measurements reveal that the asymptotic dependence of these solutions on e is that of a k n o w n class of constrained equations, then one would have evidence to suggest not only that there is in fact an underlying constrained system present (an implication which would hold even in the very weakest case: that all generic constrained systems of a given d i m e n s i o n have a c o m m o n signature), but also some of the details of its structure, such as the topology of M. Descriptive return maps could then be designed in a much less ad-hoc fashion than would otherwise be the case. This work is intended as the first step in such a classification program of constrained systems by growth law. It is an open question as to which features of the constrained system are responsible for the growth law (4) above. It may for instance be purely a function of the geometry of the constraint surface, M, rather than of any feature of the flow on it, except perhaps the n u m b e r of j u m p s trajectories undergo before returning to the section. Further work is required to answer such questions. No attempt was made to tackle this problem analytically, and it may well be possible to attain the above law by those means. For 326

13 December 1993

constrained equations of sufficient complexity however, a numerical approach such as the one presented here may be the only practical alternative. The generalization of this approach to more complex constrained equations is straightforward.

References

[ 1] B. van der Pol, Philos. Mag. 2 ( 1926) 978. [2] F. Takens, in: Lecture notes in mathematics, Vol. 525. Constrained equations, a study of implicit differential equations and their discontinuous solutions (Springer, Berlin, 1976) pp. 143-234. [3] F. Takens, in: Lecture notes in mathematics, Vol. 535. Implicit differential equations; some open problems (Springer, Berlin, 1976) pp. 237-253. [4 ] M. Diener and T. Poston, in: Springer series in synergetics. On the perfect delay convention or the revolt of the slaved variableschaos and order in nature, ed. H. Haken (Springer, Berlin, 1982). [5] O.E. Rrssler, Bull. Math. Biol. 39 (1977) 275. [6] O.E. Rrssler, Z. Naturforsch. 38a (1983) 788. [7] J.J. Tyson, J. Math. Biol. 5 (1978) 351. [8] J. Rinzel and W.C. Troy, J. Chem. Phys. 76 (1982) 1775. [9 ] J.L. Callot, F. Diener and M. Diener, C. R. Acad. Sci. A 286 (1978) 1059. [ 10] E. Benoit, J.L. Callot, F. Diener and M. Diener, Collect. Math. 32 (1981) 37. [ 11] A. Robert, Nonstandard analysis (Wiley,New York, 1988). [ 12] R. Lutz and M. Goze, in: Lecture notes in mathematics, Vol. 881. Non-standardanalysis:practicalguidewith applications (Springer, Berlin, 1981 ) pp. 399-421. [ 13] E. Benoit, Astrrisque 109/110 ( 1983) 159. [14] E. Benoit, Publ. Math. Inst. Hattes l~tudesSci. 72 (1990) 63. [ 15] F. Takens, Transitions from periodic to strange attractors in constrained equations, preprint, Groningen (1983). [ 16] E. Benoit and C. Lobry, C. R. Acad. Sci. 1294 (1982) 483. [ 17 ]L.M. Pismen, Phys. Lett. A 89 (1982) 59. [ 18 ] C. Lobryand F. Mazat, On the solutionsof some differential equationsin ~3 with smallparameter, preprint, Nice (1982). [ 19] D. Veitch, Bimodalmaps, constrainedflows, and resonance, Doctoral Thesis, University of Cambridge (1990). [20] R.S. MacKay and C. Tresser, Physica D 19 (1986 ) 206. [21]J. Milnor and W. Thurston, in: Lecture notes in mathematics, Vol. 1342. Dynamical systems, ed. J.C. Alexander (Springer, Berlin, 1988) pp. 465-563. [22] R.S. MacKay and C. Tresser, J. London Math. Soc. 37 (1988) 164. [23] G. Hall and J.M. Watt, Modem numerical methods for ordinary differentialequations (Clarendon,Oxford, 1976). [24] J. Guckenheimer and P. Holmes, Nonlinear oscillations, dynamical systems and bifurcations of vector fields (Springer, Berlin, 1983). [25 ] H. Simoyi, A. Wolf and H.L. Swinney,Phys. Rev. Lett. 49 (1982) 245.