Published in Nature, Vol. 261, p.459, June 10 1976.

Abstract. First-order difference equations arise in many contexts in the biological, economic and social sciences. Such equations, even though simple and deterministic, can exhibit a surprising array of dynamical behaviour, from stable points, to a bifurcating hierarchy of stable cycles, to apparently random fluctuations. There are consequently many fascinating problems, some concerned with delicate mathematical aspects of. the fine structure of the trajectories, and some concerned with the practical implications and applications. This is an interpretive review of them.

* King's College Research Centre, Cambridge CB2 1ST; on leave front Biology Department, Princeton University, Princeton 08540.

1. INTROUCTION

There are many situations, in many disciplines, which can be described, at least to a crude first approximation, by a simple first-order difference equation. Studies of the dynamical properties of such models usually consist of finding constant equilibrium solutions, and then conducting a linearised analysis to determine their stability with respect to small disturbances: explicitly nonlinear dynamical features are usually not considered.

Recent studies have, however, shown that the very simplest nonlinear difference equations can possess an extraordinarily rich spectrum of dynamical behaviour, from stable points, through cascades of stable cycles, to a regime in which the behaviour (although fully deterministic) is in many respects "chaotic", or indistinguishable from the sample function of a random process.

This review article has several aims.

First, although the main features of these nonlinear phenomena have been discovered and independently rediscovered by several people, I know of no source where all the main results are collected together. I have therefore tried to give such a synoptic account. This is done in a brief and descriptive way, and includes some new material: the detailed mathematical proofs are to be found in the technical literature, to which signposts are given.

Second, I indicate some of the interesting mathematical questions which do not seem to be fully resolved, Some of these problems are of a practical kind, to do with providing a probabilistic desciiption for trajectories which seem random, even though their underlying structure is deterministic. Other problems are of intrinsic mathematical interest, and treat such things as the pathology of the bifurcation structure, or the truly random behaviour, that can arise when the nonlinear function F(X) of equation (1) is not analytical. One aim here is to stimulate research on these questions, particularly on the empirical questions which relate to processing data.

Third, consideration is given to some fields where these notions may find practical application. Such applications range from the abstractly metaphorical (where, for example, the transition from a stable point to "chaos" serves as a metaphor for the onset of turbulence in a fluid), to models for the dynamic behaviour of biological populations (where one can seek to use field or laboratory data to estimate the values of the parameters in the difference equation).

Fourth, there is a very brief review of the literature pertaining to the way this spectrum of behaviour - stable points, stable cycles, chaos - can arise in second or higher order difference equations (that is, two or more dimensions; two or more interacting species), where the onset of chaos usually requires less severe nonlinearities. Differential equations are also surveyed in this light; it seems that a three-dimensional system of first-order ordinary differential equations is required for the manifestation of chaotic behaviour.

The review ends with an evangelical plea for the introduction of these difference equations into elementary mathematics courses, so that students' intuition may be enriched by seeing the wild things that simple nonlinear equations can do.

2. FIRST-ORDER DIFFERENCE EQUATIONS

One of the simplest systems an ecologist can study is a seasonally breeding population in which generations do not overlap ^{1 - 4}. Many natural populations, particularly among temperate zone insects (including many economically important crop and orchard pests), are of this kind, In this situation, the observational data will usually consist of information about the maximum, or the average, or the total population in each generation. The theoretician seeks to understand how the magnitude of the population in generation t+1, X_t+1, is related to the magnitude of the population in the preceding generation t, X_t: such a relationship may be expressed in the general form

The function F(X) will usually be what a biologist calls "density dependent", and a mathematician calls nonlinear; equation (1) is then a first-order, nonlinear difference equation.

Although I shall henceforth adopt the habit of referring to the variable X as "the population", there are countless situations outside population biology where the basic equation (1), applies. There are other examples in biology, as, for example in genetics ^{5, 6} (where the equation describes the change in gene frequency in time) or in epidemiology ⁷ (with X the fraction of the population infected at time t). Examples in economics include models for the relationship between commodity quantity and price ⁸, for the theory of business cycles ⁹, and for the temporal sequences generated by various other economic quantities ¹⁰. The general equation (1) also is germane to the social sciences ¹¹, where it arises, for example, in theories of learning (where X may be the number of bits of information that can be remembered after an interval t), or in the propagation of rumours in variously structured societies (where X is the number of people to have heard the rumour after time t). The imaginative reader will be able to invent other contexts for equation (1).

In many of these contexts, and for biological populations in particular, there is a tendency for the variable X to increase from one generation to the next when it is small, and for it to decrease when it is large. That is, the nonlinear function F(X) often has the following properties: F(0)=0; F(X) increases monotonically as X increases through the range 0 < X < A (with F(X) attaining its maximum value at X = A); and F(X) decreases monotonically as X increases beyond X = A. Moreover, F(X) will usually contain one or more parameters which "tune" the severity of this nonlinear behaviour; parameters which tune the steepness of the hump in the F(X) curve. These parameters will typically have some biological or economic or sociological significance.

A specific example is afforded by the equation ^{1, 4, 12 - 23}

This is sometimes called the "logistic" difference equation. In the limit b = 0, it describes a population growing purely exponentially (for a > 1); for b 0, the quadratic nonlinearity produces a growth curve with a hump, the steepness of which is tuned by the parameter a. By writing X = bN/a, the equation may be brought into canonical form ^{1, 4, 12 - 23}

In this form, which is illustrated in Fig. 1, it is arguably the simplest nonlinear difference equation. I shall use equation (3) for most of the numerical examples and illustrations in this article. Although attractive to mathematicians by virtue of its extreme simplicity, in practical applications equation (3) has the disadvantage that it requires X to remain on the interval 0 < X < 1; if X ever exceeds unity, subsequent iterations diverge towards - (which means the population becomes extinct). Furthermore, F(X) in equation (3) attains a maximum value of a/4 (at X = 1/2); the equation therefore possesses non-trivial dynamical behaviour only if a < 4. On the other hand, all trajectories are attracted to X = 0 if a < 1. Thus for non-trivial dynamical behaviour we require 1 < a < 4; failing this, the population becomes extinct.

Figure 1. A typical form for the relationship between Xt+1 and Xt described by equation (1). The curves are for equation (3), with a = 2.707 (a); and a = 3.414 (b). The dashed lines indicate the slope at the "fixed points" where F(X) intersects the 45° line: for the case a this slope is less steep than -45° and the fixed point is stable; for b the slope is steeper than -45°, and the point is unstable.

Another example, with a more secure provenance in the biological literature 1, 23 - 27, is the equation

This again describes a population with a propensity to simple exponential growth at low densities, and a tendency to decrease at high densities. The steepness of this nonlinear behaviour is tuned by the parameter r. The model is plausible for a single species population which is regulated by an epidemic disease at high density ²⁸. The function F(X) of equation (4) is slightly more complicated than that of equation (3), but has the compensating advantage that local stability implies global stability¹ for all X > 0.

The forms (3) and (4) by no means exhaust the list of single-humped functions F(X) for equation (1) which can be culled from the ecological literature. A fairly full such catalogue is given, complete with references, by May and Oster ¹. Other similar mathematical functions are given by Metropolis et al. ¹⁶. Yet other forms for F(X) are discussed under the heading of "mathematical curiosities" below.

3. DYNAMIC PROPERTIES OF EQUATION (1)

Possible constant, equilibrium values (or "fixed points") of X in equation (1) may be found algebraically by putting X_t+1 = X_t = X*, and solving the resulting equation

An equivalent graphical method is to find the points where the curve F(X) that maps X_t into X_t+1 intersects the 45° line, X_t+1 = X_t which corresponds to the ideal nirvana of zero population growth; see Fig. 1. For the single-hump curves discussed above, and exemplified by equations (3) and (4), there are two such points: the trivial solution X = 0, and a non-trivial solution X* (which for equation (3) is X* = 1 - (1/a).

The next question concerns the stability of the equilibrium point X*. This can be seen ^{24, 25, 19 - 21, 1, 4} to depend on the slope of the F(X) curve at X*. This slope, which is illustrated by the dashed lines in Fig. 1, can be designated

(X*) = [dF / dX]_{x = x*}         (6)

So long as this slope lies between 45° and -45° (that is, ⁽¹⁾ between +1 and -1), making an acute angle with the 45° ZPG line, the equilibrium point X* will be at least locally stable, attracting all trajectories in its neighbourhood. In equation (3), for example, this slope is ⁽¹⁾ = 2 - a: the equilibrium point is therefore stable, and attracts all trajectories originating in the interval 0 < X < 1, if and only if 1 < a < 3.

As the relevant parameters are tuned so that the curve F(X) becomes more and more steeply humped, this stability-determining slope at X* may eventually steepen beyond -45° (that is, ⁽¹⁾ < -1), whereupon the equilibrium point X* is no longer stable.

What happens next? What happens, for example, for a > 3 in equation (3)?

To answer this question, it is helpful to look at the map which relates the populations at successive intervals 2 generations apart; that is, to look at the function which relates X_t+2 to X_t. This second iterate of equation (1) can be written

or, introducing an obvious piece of notation,

The map so derived from equation (3) is illustrated in Figs 2 and 3.

Figure 2. The map relating Xt+2 to Xt, obtained by two iterations of equation (3). This figure is for the case (a) of Fig. 1, a = 2.707: the basic fixed point is stable, and it is the only point at which F(2)(X) intersects the 45° line (where its slope, shown by the dashed line, is less steep than 45°).

Figure 3. As for Fig. 2, except that here a = 3.414, as in Fig. 1b. The basic fixed point is now unstable: the slope of F(2)(X) at this point steepens beyond 45°, leading to the appearance of two new solutions of period 2.

or graphically from the intersection between the map F⁽²⁾(X) and the 45° line, as shown in Figs 2 and 3. Clearly the equilibrium point X* of equation (5) is a solution of equation (9); the basic fixed point of period 1 is a degenerate case of a period 2 solution. We now make a simple, but crucial, observation ¹: the slope of the curve F⁽²⁾(X) at the point X*, defined as ⁽²⁾(X*) and illustrated by the dashed lines in Figs 2 and 3, is the square of the corresponding slope of F(X)

⁽²⁾ (X*) = [⁽¹⁾ (X*)]²         (10)

This fact can now be used to make plain what happens when the fixed point X* becomes unstable. If the slope of F(X) is less than -45° (that is, |⁽¹⁾| < 1), as illustrated by curve a in Fig. 1, then X* is stable. Also, from equation (10), this implies 0 < ⁽²⁾ < 1 corresponding to the slope of F⁽²⁾ at X* lying between 0° and 45°, as shown in Fig. 2. As long as the fixed point X* is stable, it provides the only non-trivial solution to equation (9). On the other hand, when ⁽¹⁾ steepens beyond -45° (that is, |⁽¹⁾| > 1), as illustrated by curve b in Fig 1, X* becomes unstable. At the same time, from equation (10) this implies ⁽²⁾ > 1, corresponding to the slope of F⁽²⁾ at X* steepening beyond 45°, as shown in Fig. 3. As this happens, the curve F⁽²⁾(X) must develop a "loop", and two new fixed points of period 2 appear, as illustrated in Fig. 3.

In short, as the nonlinear function F(X) in equation (1) becomes more steeply humped, the basic fixed point X* may become unstable. At exactly the stage when this occurs, there are born two new and initially stable fixed points of period 2, between which the system alternates in a stable cycle of period 2. The sort of graphical analysis indicated by Figs 1, 2 and 3, along with the equation (10), is all that is needed to establish this generic result ^{1, 4}.

As before, the stability of this period 2 cycle depends on the slope of the curve F⁽²⁾(X) at the 2 points. (This slope is easily shown to be the same at both points ^{1, 20}, and more generally to be the same at all k points on a period k cycle.) Furthermore, as is clear by imagining the intermediate stages between Figs 2 and 3, this stability-determining slope has the value = +1 at the birth of the 2-point cycle, and then decreases through zero towards = -1 as the hump in F(X) continues to steepen. Beyond this point the period 2 points will in turn become unstable, and bifurcate to give an initially stable cycle of period 4. This in turn gives way to a cycle of period 8, and thence to a hierarchy of bifurcating stable cycles of periods 16, 32, 64,..., 2ⁿ. In each case, the way in which a stable cycle of period k becomes unstable, simultaneously bifurcating to produce a new and initially stable cycle of period 2k, is basically similar to the process just adumbrated for k = 1. A more full and rigorous account of the material covered so far is in ref. 1.

This "very beautiful bifurcation phenomenon" ²² is depicted in Fig. 4, for the example equation (3). It cannot be too strongly emphasised that the process is generic to most functions F(X) with a hump of tunable steepness. Metropolis et al. ¹⁶ refer to this hierarchy of cycles of periods 2ⁿ as the harmonics of the fixed point X*.

Figure 4. This figure illustrates some of the stable ( ___ ) and unstable (----) fixed points of various periods that can arise by bifurcation processes in equation (1) in general, and equation (3) in particular. To the left, the basic stable fixed point becomes unstable and gives rise by a succession of pitchfork bifurcations to stable harmonics of period 2n; none of these cycles is stable beyond a = 3.5700. To the right, the two period 3 cycles appear by tangent bifurcation: one is initially unstable; the other is initially stable, but becomes unstable and gives way to stable harmonics of period 3 × 2n, which have a point of accumulation at a = 3.8495. Note the change in scale on the a axis, needed to put both examples on the same figure. There are infinitely many other such windows, based on cycles of higher periods.

Beyond this point of accumulation (for example, for a > a_c in equation (3)) there are an infinite number of fixed points with different periodicities, and an infinite number of different periodic cycles. There are also an uncountable number of initial points X₀ which give totally aperiodic (although bounded) trajectories; no matter how long the time series generated by F(X) is run out, the pattern never repeats. These facts may be established by a variety of methods 1, 4, 20, 22, 29. Such a situation, where an infinite number of different orbits can occur, has been christened "chaotic" by Li and Yorke²⁰.

As the parameter increases beyond the critical value, at first all these cycles have even periods, with X_t alternating up and down between values above, and values below, the fixed point X*. Although these cycles may in fact be very complicated (having a non-degenerate period of, say, 5,726 points before repeating), they will seem to the casual observer to be rather like a somewhat "noisy" cycle of period 2. As the parameter value continues to increase, there comes a stage (at a = 3.6786.. for equation (3)) at which the first odd period cycle appears. At first these odd cycles have very long periods, but as the parameter value continues to increase cycles with smaller and smaller odd periods are picked up, until at last the three-point cycle appears (at a = 3.8284 . . for equation (3)). Beyond this point, there are cycles with every integer period, as well as an uncountable number of asymptotically aperiodic trajectories: Li and Yorke ²⁰ entitle their original proof of this result "Period Three Implies Chaos".

The term "chaos" evokes an image of dynamical trajectories which are indistinguishable from some stochastic process. Numerical simulations ^{12, 15, 21, 23, 25} of the dynamics of equation (3), (4) and other similar equations tend to confirm this impression. But, for smooth and "sensible" functions F(X) such as in equations (3) and (4), the underlying mathematical fact is that for any specified parameter value there is one unique cycle that is stable, and that attracts essentially all initial points ^{22, 29} (see ref. 4, appendix A, for a simple and lucid exposition). That is, there is one cycle that "owns" almost all initial points; the remaining infinite number of other cycles, along with the asymptotically aperiodic trajectories, own a set of points which, although uncountable, have measure zero.

As is made clear by Tables 3 and 4 below, any one particular stable cycle is likely to occupy an extraordinarily narrow window of parameter values. This fact, coupled with the long time it is likely to take for transients associated with the initial conditions to damp out, means that in practice the unique cycle is unlikely to be unmasked, and that a stochastic description of the dynamics is likely to be appropriate, in spite of the underlying deterministic structure. This point is pursued further under the heading "practical applications", below.

The main messages of this section are summarised in Table 1, which sets out the various domains of dynamical behaviour of the equations (3) and (4) as functions of the parameters, a and r respectively, that determine the severity of the nonlinear response. These properties can be understood qualitatively in a graphical way, and are generic to any well behaved F(X) in equation (1).

Table 1. Summary of the way various "single-hump" functions F(X), from equation (1), behave in the chaotic region, distinguishing the dynamical properties which are generic from those which are not

The function F(X)			aX; if X < ½	X; if X < 1
of equation (1)	aX(1 - X)	X exp[r(1 - X)]	a(1 - X); if X > ½	X1-b; if X > 1
Tunable parameter	a	r	a	b
Fixed point becomes unstable	3.0000	2.0000	1.0000*	2.0000
"Chaotic" region begins
[point of accumulation of cycles of period 2n]	3.5700	2.6924	1.0000	2.0000
First odd-period cycle appears	3.6786	2.8332	1.4142	2.6180
Cycle with period 3 appears
[and therefore every integer period present]	3.8284	3.1024	1.6180	3.0000
"Chaotic" region ends	4.0000†	‡	2.000†	‡
Are there stable cycles in the chaotic region?	Yes	Yes	No	No
* Below this a value, X = 0 is stable.
† All solutions are attracted to - for a values beyond this.
‡ In practice, as r or b becomes large enough, X will eventually be carried so low as to be effectively zero, thus producing extinction in models of biological populations.

We now proceed to a more detailed discussion of the mathematical structure of the chaotic regime for analytical functions, and then to the practical problems alluded to above and to a consideration of the behavioural peculiarities exhibited by non-analytical functions (such as those in the two right hand columns of Table 1).

4. FINE STRUCTURE OF THE CHAOTIC REGIME

We have seen how the original fixed point X* bifurcates to give harmonics of period 2ⁿ. But how do new cycles of period k arise?

The general process is illustrated in Fig. 5, which shows how period 3 cycles originate. By an obvious extension of the notation introduced in equation (8). populations three generations apart are related by

If the hump in F(X) is sufficiently steep, the threefold iteration will produce a function F⁽³⁾(X) with 4 humps, as shown in Fig. 5 for the F(X) of equation (3). At first (for a < 3.8284 . . in equation 3) the 45° line intersects this curve only at the single point X* (and at X = 0), as shown by the solid curve in Fig. 5. As the hump in F(X) steepens, the hills and valleys in F⁽³⁾(X) become more pronounced, until simultaneously the first two valleys sink and the final hill rises to touch the 45° line, and then to intercept it at 6 new points, as shown by the dashed curve in Fig. 5. These 6 points divide into two distinct three-point cycles. As can be made plausible by imagining the intermediate stages in Fig. 5, it can be shown that the stability-determining slope of F⁽³⁾(X) at three of these points has a common value, which is ⁽³⁾ = +1 at their birth, and thereafter steepens beyond +1: this period 3 cycle is never stable. The slope of F⁽³⁾(X) at the other three points begins at ⁽³⁾ = +1, and then decreases towards zero, resulting in a stable cycle of period 3. As F(X) continues to steepen, the slope ⁽³⁾ for this initially stable three-point cycle decreases beyond -1; the cycle becomes unstable, and gives rise by the bifurcation process discussed in the previous section to stable cycles of period 6, 12, 24, ..., 3 × 2ⁿ. This birth of a stable and unstable pair of period 3 cycles, and the subsequent harmonics which arise as the initially stable cycle becomes unstable, are illustrated to the right of Fig. 4.

Figure 5. The relationship between Xt+3 and Xt, obtained by three iterations of equation (3). The solid curve is for a = 3.7, and only intersects the 45° line once. As a increases, the hills and valleys become more pronounced. The dashed curve is for a = 3.9, and six new period 3 points have appeared (arranged as two cycles, each of period 3).

There are, therefore, two basic kinds of bifurcation processes 1, 4 for first order difference equations. Truly new cycles of period k arise in pairs (one stable, one unstable) as the hills and valleys of higher iterates of F(X) move, respectively, up and down to intercept the 45° line, as typified by Fig. 5. Such cycles are born at the moment when the hills and valleys become tangent to the 45° line, and the initial slope of the curve F^(k) at the points is thus ^(k) = +1: this type of bifurcation may be called ^{1, 4} a tangent bifurcation or a = +1 bifurcation. Conversely, an originally stable cycle of period k may become unstable as F(X) steepens. This happens when the slope of F^(k) at these period k points steepens beyond ^(k) = -1, whereupon a new and initially stable cycle of period 2k is born in the way typified by Figs 2 and 3. This type of bifurcation may be called a pitchfork bifurcation (borrowing an image from the left hand side of Fig. 4) or a = -1 bifurcation ^{1, 4}.

Putting all this together, we conclude that as the parameters in F(X) are varied the fundamental, stable dynamical units are cycles of basic period k, which arise by tangent bifurcation, along with their associated cascade of harmonics of periods k2ⁿ, which arise by pitchfork bifurcation. On this basis, the constant equilibrium solution X* and the subsequent hierarchy of stable cycles of periods 2ⁿ is merely a special case, albeit a conspicuously important one (namely k = 1), of a general phenomenon. In addition, remember ^{1, 4, 22, 29} that for sensible, analytical functions (such as, for example, those in equations (3) and (4)) there is a unique stable cycle for each value of the parameter in F(X). The entire range of parameter values (1 < a < 4 in equation (3), 0 < r in equation (4)) may thus be regarded as made up of infinitely many windows of parameter values - some large, some unimaginably small - each corresponding to a single one of these basic dynamical units. Tables 3 and 4, below, illustrate this notion. These windows are divided from each other by points (the points of accumulation of the harmonics of period k2ⁿ) at which the system is truly chaotic, with no attractive cycle: although there are infinitely many such special parameter values, they have measure zero on the interval of all values.

How are these various cycles arranged along the interval of relevant parameter values? This question has to my knowledge been answered independently by at least 6 groups of people, who have seen the problem in the context of combinatorial theory ^{16, 30}, numerical analysis^{13, 14}, population biology ¹, and dynamical systems theory ^{22, 31} (broadly defined).

A simple-minded approach (which has the advantage of requiring little technical apparatus, and the disadvantage of being rather clumsy) consists of first answering the question, how many period k points can there be? That is, how many distinct solutions can there be to the equation

If the function F(X) is sufficiently steeply humped, as it will be once the parameter values are sufficiently large, each successive iteration doubles the number of humps, so that F^(k)(X) has 2^k-1 humps. For large enough parameter values, all these hills and valleys will intersect the 45° line, producing 2^k fixed points of period k. These are listed for k 12 in the top row of Table 2. Such a list includes degenerate points of period k, whose period is a submultiple of k; in particular, the two period 1 points (X = 0 and X*) are degenerate solutions of equation (12) for all k. By working from left to right across Table 2, these degenerate points can be subtracted out, to leave the total number of non-degenerate points of basic period k, as listed in the second row of Table 2. More sophisticated ways of arriving at this result are given elsewhere ^{13, 14, 16, 22, 30, 31}.

The 2^k period k points are arranged into various cycles of period k, or submultiples thereof, which appear in succession by either tangent or pitchfork bifurcation as the parameters in F(X) are varied. The third row in Table 2 catalogues the total number of distinct cycles of period k which so appear. In the fourth row ¹⁴, the degenerate cycles are subtracted out, to give the total number of non-degenerate cycles of period k: these numbers must equal those of the second row divided by k. This fourth row includes the (stable) harmonics which arise by pitchfork bifurcation, and the pairs of stable-unstable cycles arising by tangent bifurcation. By subtracting out the cycles which are unstable from birth, the total number of possible stable cycles is given in row five; these figures can also be obtained by less pedestrian methods ^{13, 16, 30}. Finally we may subtract out the stable cycles which arise by pitchfork bifurcation, as harmonics of some simpler cycle, to arrive at the final row in Table 2, which lists the number of stable cycles whose basic period is k.

Returning to the example of period 6, we have already noted the five degenerate cycles whose periods are submultiples of 6. The remaining 54 points are parcelled out into one cycle of period 6 which arises as the harmonic of the only stable three-point cycle, and four distinct pairs of period 6 cycles (that is, four initially stable ones and four unstable ones) which arise by successive tangent bifurcations. Thus, reading from the foot of the column for period 6 in Table 2, we get the numbers 4, 5, 9, 14.

Using various labelling tricks, or techniques from combinatorial theory, it is also possible to give a generic list of the order in which the various cycles appear ^{1, 13, 16, 22}. For example, the basic stable cycles of periods 3, 5, 6 (of which there are respectively 1, 3, 4) must appear in the order 6, 5, 3, 5, 6, 6, 5, 6: compare Tables 3 and 4. Metropolis et al. ¹⁶ give the explicit such generic list for all cycles of period k 11.

As a corollary it follows that, given the most recent cycle to appear, it is possible (at least in principle) to catalogue all the cycles which have appeared up to this point. An especially elegant way of doing this is given by Smale and Williams 22, who show, for example, that when the stable cycle of period 3 first originates, the total number of other points with periods k, N_k, which have appeared by this stage satisfy the Fibonacci series, N_k = 2, 4, 5, 8, 12, 19, 30, 48, 77, 124, 200, 323 for k = 1, 2, ..., 12: this is to be contrasted with the total number of points of period k which will eventually appear (the top row of Table 2) as F(X) continues to steepen.

Hoppensteadt and Hyman 21 have given a simple graphical method for locating the parameter value in the chaotic regime at which the first odd period cycle appears. Their analytic recipe is as follows. Let be the parameter which tunes the steepness of F(X) (for example, = a for equation (3), = r for equation (4)), X*() be the fixed point of period 1 (the non-trivial solution of equation (5)), and X_max() the maximum value attainable from iterations of equation (1) (that is, the value of F(X) at its hump or stationary point). The first odd period cycle appears for that value of which satisfies ^{21, 31}

X*() = F⁽²⁾ (X_max ())         (13)

As mentioned above, another critical value is that where the period 3 cycle first appears. This parameter value may be found numerically from the solutions of the third iterate of equation (1): for equation (3) it is ¹⁴ a = 1 + 8.

Myrberg ¹³ (for all k 10) and Metropolis et al. ¹⁶. (for all k 7) have given numerical information about the stable cycles in equation (3). They do not give the windows of parameter values, but only the single value at which a given cycle is maximally stable; that is, the value of a for which the stability-determining slope of F^(k)(X) is zero, ^(k) = 0. Since the slope of the k-times iterated map F^(k) at any point on a period k cycle is simply equal to the product of the slopes of F(X) at each of the points X*_k on this cycle ^{1, 8, 20}, the requirement ^(k) = 0 implies that X = A (the stationary point of F(X), where ⁽¹⁾ = 0) is one of the periodic points in question, which considerably simplifies the numerical calculations.

For each basic cycle of period k (as catalogued in the last row of Table 2), it is more interesting to know the parameter values at which: (1) the cycle first appears (by tangent bifurcation); (2) the basic cycle becomes unstable (giving rise by successive pitchfork bifurcations to a cascade of harmonics of periods k2ⁿ); (3) all the harmonics become unstable (the point of accumulation of the period k2ⁿ cycles). Tables 3 and 4 extend the work of May and Oster¹, to give this numerical information for equations (3) and (4), respectively. (The points of accumulation are not ground out mindlessly, but are calculated by a rapidly convergent iterative procedure, see ref. 1, appendix A.) Some of these results have also been obtained by Gumowski and Mira ³².

5. PRACTICAL PROBLEMS

Referring to the paradigmatic example of equation (3), we can now see that the parameter interval 1 < a < 4 is made up of a one-dimensional mosaic of infinitely many windows of a-values, in each of which a unique cycle of period k, or one of its harmonics, attracts essentially all initial points. Of these windows, that for 1 < a < 3.5700 . . corresponding to k = 1 and its harmonics is by far the widest and most conspicuous. Beyond the first point of accumulation, it can be seen from Table 3 that these windows are narrow, even for cycles of quite low periods, and the windows rapidly become very tiny as k increases.

As a result, there develops a dichotomy between the underlying mathematical behaviour (which is exactly determinable) and the "commonsense" conclusions that one would draw from numerical simulations. If the parameter a is held constant at one value in the chaotic region, and equation (3) iterated for an arbitrarily large number of generations, a density plot of the observed values of X_t on the interval 0 to 1 will settle into k equal spikes (more precisely, delta functions) corresponding to the k points on the stable cycle appropriate to this a-value. But for most a-values this cycle will have a fairly large period, and moreover it will typically take many thousands of generations before the transients associated with the initial conditions are damped out: thus the density plot produced by numerical simulations usually looks like a sample of points taken from some continuous distribution.

An especially interesting set of numerical computations are due to Hoppensteadt (personal communication) who has combined many iterations to produce a density plot of X_t for each one of a sequence of a-values, gradually increasing from 3.5700 . . to 4. These results are displayed as a movie. As can be expected from Table 3, some of the more conspicuous cycles do show up as sets of delta functions: the 3-cycle and its first few harmonics; the first 5-cycle; the first 6-cycle. But for most values of a the density plot looks like the sample function of a random process. This is particularly true in the neighbourhood of the a-value where the first odd cycle appears (a = 3.6786..), and again in the neighbourhood of a = 4: this is not surprising, because each of these locations is a point of accumulation of points of accumulation. Despite the underlying discontinuous changes in the periodicities of the stable cycles, the observed density pattern tends to vary smoothly. For example, as a increases toward the value at which the 3-cycle appears, the density plot tends to concentrate around three points, and it smoothly diffuses away from these three points after the 3-cycle and all its harmonics become unstable.

I think the most interesting mathematical problem lies in designing a way to construct some approximate and "effectively continuous" density spectrum, despite the fact that the exact density function is determinable and is always a set of delta functions. Perhaps such techniques have already been developed in ergodic theory ³³ (which lies at the foundations of statistical mechanics), as for example in the use of "coarse-grained observers". I do not know.

Such an effectively stochastic description of the dynamical properties of equation (4) for large r has been provided ²⁸, albeit by tactical tricks peculiar to that equation rather than by any general method. As r increases beyond about 3, the trajectories generated by this equation are, to an increasingly good approximation, almost periodic with period (1 / r) exp(r - 1).

The opinion I am airing in this section is that although the exquisite fine structure of the chaotic regime is mathematically fascinating, it is irrelevant for most practical purposes. What seems called for is some effectively stochastic description of the deterministic dynamics. Whereas the various statements about the different cycles and their order of appearance can be made in generic fashion, such stochastic description of the actual dynamics will be quite different for different F(X): witness the difference between the behaviour of equation (4), which for large r is almost periodic "outbreaks" spaced many generations apart, versus the behaviour of equation (3), which for a -> 4 is not very different from a series of Bernoulli coin flips.

6. MATHEMATICAL CURIOSITIES

As discussed above, the essential reason for the existence of a succession of stable cycles throughout the "chaotic" regime is that as each new pair of cycles is born by tangent bifurcation (see Fig. 5), one of them is at first stable, by virtue of the way the smoothly rounded hills and valleys intercept the 45° line. For analytical functions F(X), the only parameter values for which the density plot or "invariant measure" is continuous and truly ergodic are at the points of accumulation of harmonics, which divide one stable cycle from the next. Such exceptional parameter values have found applications, for example, in the use of equation (3) with a = 4 as a random number generator ^{34, 35}: it has a continuous density function proportional to [X(1 - X)]^-1/2 in the interval 0 < X < 1.

Non-analytical functions F(X) in which the hump is in fact a spike provide an interesting special case. Here we may imagine spikey hills and valleys moving to intercept the 45° line in Fig. 5, and it may be that both the cycles born by tangent bifurcation are unstable from the outset (one having ^(k) > 1, the other ^(k) < -1), for all k > 1. There are then no stable cycles in the chaotic regime, which is therefore literally chaotic with a continuous and truly ergodic density distribution function.

One simple example is provided by

defined on the interval 0 < X < 1. For 0 < a < 1, all trajectories are attracted to X = 0; for 1 < a < 2, there are infinitely many periodic orbits, along with an uncountable number of aperiodic trajectories, none of which are locally stable. The first odd period cycle appears at a = 2, and all integer periods are represented beyond a = (1 + 5)/2. Kac³⁶ has given a careful discussion of the case a = 2. Another example, this time with an extensive biological pedigree ^{1 - 3}, is the equation

Xt+1 = Xt ; if Xt < 1

Xt+1 = Xt1-b ; if Xt > 1         (15)

If > 1 this possesses a globally stable equilibrium point for b < 2. For b > 2 there is again true chaos, with no stable cycles: the first odd cycle appears at b = (3 + 5)/2, and all integer periods are present beyond b = 3. The dynamical properties of equations (14) and (15) are summarised to the right of Table 2.

The absence of analyticity is a necessary, but not a sufficient, condition for truly random behaviour ³¹. Consider, for example,

This is the parabola of equation (3) and Fig. 1, but with the left hand half of F(X) flattened into a straight line. This equation does possess windows of a values, each with its own stable cycle, as described generically above. The stability-determining slopes ^(k) vary, however, discontinuously with the parameter a, and the widths of the simpler stable regions are narrower than for equation (3): the fixed point becomes unstable at a = 3; the point of accumulation of the subsequent harmonics is at a = 3.27 . .; the first odd cycle appears at a = 3.44 . .; the 3-point cycle at a = 3.67. . (compare the first column in Table 1).

These eccentricities of behaviour manifested by non-analytical functions may be of interest for exploring formal questions in ergodic theory. I think, however, that they have no relevance to models in the biological and social sciences, where functions such as F(X) should be analytical. This view is elaborated elsewhere ³⁷.

As a final curiosity, consider the equation

X_t+1 = X_t [1 + X_t]^-         (17)

This has been used to fit a considerable amount of data on insect populations ^{38, 39}. Its stability behaviour, as a function of the two parameters and , is illustrated in Fig. 6. Notice that for < 7.39 . . there is a globally stable equilibrium point for all ; for 7.39. . < X < 12.50. . this fixed point becomes unstable for sufficiently large , bifurcating to a stable 2-point cycle which is the solution for all larger ; as increases through the range 12.50 . . < < 14.77 . . various other harmonics of period 2ⁿ appear in turn. The hierarchy of bifurcating cycles of period 2ⁿ is thus truncated, and the point of accumulation and subsequent regime of chaos is not achieved (even for arbitrarily large ) until > 14.77...

Figure 6. The solid lines demarcate the stability domains for the density dependence parameter, , and the population growth rate, , in equation (17); the dashed line shows where 2-point cycles give way to higher cycles of period 2n. The solid circles come from analyses of life table data on field populations, and the open circles from laboratory populations (from ref. 3, after ref. 39).

The fact that the simple and deterministic equation (1) can possess dynamical trajectories which look like some sort of random noise has disturbing practical implications. It means, for example, that apparently erratic fluctuations in the census data for an animal population need not necessarily betoken either the vagaries of an unpredictable environment or sampling errors: they may simply derive from a rigidly deterministic population growth relationship such as equation (1). This point is discussed more fully and carefully elsewhere 1.

Alternatively, it may be observed that in the chaotic regime arbitrarily close initial conditions can lead to trajectories which, after a sufficiently long time, diverge widely. This means that, even if we have a simple model in which all the parameters are determined exactly, long term prediction is nevertheless impossible. In a meteorological context, Lorenz ¹⁵ has called this general phenomenon the "butterfly effect": even if the atmosphere could be described by a deterministic model in which all parameters were known, the fluttering of a butterfly's wings could alter the initial conditions, and thus (in the chaotic regime) alter the long term prediction.

Fluid turbulence provides a classic example where, as a parameter (the Reynolds number) is tuned in a set of deterministic equations (the Navier-Stokes equations), the motion can undergo an abrupt transition from some stable configuration (for example, laminar flow) into an apparently stochastic, chaotic regime. Various models, based on the Navier-Stokes differential equations, have been proposed as mathematical metaphors for this process ^{15, 40, 41}. In a recent review of the theory of turbulence, Martin ⁴² has observed that the one-dimensional difference equation (1) may be useful in this context. Compared with the earlier models ^{15, 40, 41}, it has the disadvantage of being even more abstractly metaphorical, and the advantage of having a spectrum of dynamical behaviour which is more richly complicated yet more amenable to analytical investigation.

A more down-to-earth application is possible in the use of equation (1) to fit data ^{1, 2, 3, 38, 39, 43} on biological populations with discrete, non-overlapping generations, as is the case for many temperate zone arthropods. Figure 6 shows the parameter values and that are estimated ³⁹ for 24 natural populations and 4 laboratory populations when equation (17) is fitted to the available data. The figure also shows the theoretical stability domains: a stable point; its stable harmonics (stable cycles of period 2ⁿ); chaos. The natural populations tend to have stable equilibrium point behaviour. The laboratory populations tend to show oscillatory or chaotic behaviour; their behaviour may be exaggeratedly nonlinear because of the absence, in a laboratory setting, of many natural mortality factors. It is perhaps suggestive that the most oscillatory natural population (labelled A in Fig. 6) is the Colorado potato beetle, whose present relationship with its host plant lacks an evolutionary pedigree. These remarks are only tentative, and must be treated with caution for several reasons. Two of the main caveats are that there are technical difficulties in selecting and reducing the data, and that there are no single species populations in the natural world: to obtain a one-dimensional difference equation by replacing a population's interactions with its biological and physical environment by passive parameters (such as and ) may do great violence to the reality.

Some of the many other areas where these ideas have found applications were alluded to in the second section, above ^{5 - 11}. One aim of this review article is to provoke applications in yet other fields.

8. RELATED PHENOMENA IN HIGHER DIMENSIONS

Pairs of coupled, first-order difference equations (equivalent to a single second-order equation) have been investigated in several contexts ^{4, 44 - 46}, particularly in the study of temperate zone arthropod prey-predator systems ^{2 - 4, 23, 47}. In these two-dimensional systems, the complications in the dynamical behaviour are further compounded by such facts as: (1) even for analytical functions, there can be truly chaotic behaviour (as for equations (14) and (15)), corresponding to so-called "strange attractors"; and (2) two or more different stable states (for example, a stable point and a stable cycle of period 3) can occur together for the same parameter values ⁴. In addition, the manifestation of these phenomena usually requires less severe nonlinearities (less steeply humped F(X)) than for the one-dimensional case.

Similar systems of first-order ordinary differential equations, or two coupled first-order differential equations, have much simpler dynamical behaviour, made up of stable and unstable points and limit cycles ⁴⁸. This is basically because in continuous two-dimensional systems the inside and outside of closed curves can be distinguished; dynamic trajectories cannot cross each other. The situation becomes qualitatively more complicated and in many ways analogous to first-order difference equations when one moves to systems of three or more coupled, first-order ordinary differential equations (that is, three-dimensional systems of ordinary differential equations). Scanlon (personal communication) has argued that chaotic behaviour and "strange attractors", that is solutions which are neither points nor periodic orbits ⁴⁸, are typical of such systems. Some well studied examples arise in models for reaction-diffusion systems in chemistry and biology ⁴⁹, and in the models of Lorenz ¹⁵ (three dimensions) and Ruelle and Takens ⁴⁰ (four dimensions) referred to above. The analysis of these systems is, by virtue of their higher dimensionality, much less transparent than for equation (1).

An explicit and rather surprising example of a system which has recently been studied from this viewpoint is the ordinary differential equations used in ecology to describe competing species. For one or two species these systems are very tame: dynamic trajectories will converge on some stable equilibrium point (which may represent coexistence, or one or both species becoming extinct). As Smale ⁵⁰ has recently shown, however, for 3 or more species these general equations can, in a certain reasonable and well-defined sense, be compatible with any dynamical behaviour. Smale's ⁵⁰ discussion is generic and abstract: a specific study of the very peculiar dynamics which can be exhibited by the familiar Lotka-Volterra equations once there are 3 competitors is given by May and Leonard ⁵¹.

REFERENCES

1. May, R. M., and Oster, G. F., Am. Nat., 110 (in the press).
2. Varley, G. C., Gradwell, G. R., and Hassell, M. P., Insect Population Ecology (Blackwell, Oxford. 1973).
3. May, R. M. (ed.), Theoreticol Ecology: Principles and Applications (Blackwell, Oxford, 1976).
4. Guckenheimer, J., Oster, G. F., and Ipaktchi, A., Theor. Pop. Biol. (in the press).
5. Oster, G. F., Ipaktchi, A., and Rocklin, I., Theor. Pop. Biol. (in the press).
6. Asmussen, M. A., and Feldman, M. W., J. theor. Biol. (in the press).
7. Hoppensteadt, F. C; Mathematical Theories of Populations: Demographics, Genetics and Epidemics (SIAM, Philadelphia, 1975).
8. Samuelson, P. A., Foundations of Economic Analysis (Harvard University Press, Cambridge, Massachusetts, 1947).
9. Goodwin, R. E., Econometrica, 19, 1-17 (1951).
10. Baumol, W. J., Economnic Dynamics, 3rd ed. (Macmillan, New York, 1970).
11. See, for example, Kemeny, J., and Snell, J. L., Mathematical Models in the Social Sciences (MIT Press, Cambridge, Massachusetts, 1972).
12. Chaundy, T. W., and Phillips, E., Q. JI Math. Oxford, 7, 74-80 (1936).
13. Myrberg, P. J., Ann. Akad. Sc. Fennicae, A, I, No. 336/3 (1963).
14. Myrberg, P. J., Ann. Akad. Sc. Fennicae, A, I, No. 259 (1958).
15. Lorenz, E. N., J. Atmos. Sci., 20, 130-141 (1963); Tellus, 16, 1-11 (1964).
16. Metropolis, N., Stein, M. L., and Stein, P. R., J. Combinatorial Theory, 15(A), 25-44 (1973).
17. Maynard Smith, J., Mathematical Ideas in Biology (Cambridge University Press, Cambridge, 1968).
18. Krebs, C. J., Ecology (Harper and Row, New York, 1972).
19. May, R. M., Am. Nat., 107, 46-57 (1972).
20. Li, T-Y., and Yorke, J. A., Am. Math. Monthly, 82, 985-992 (1975).
21. Hoppensteadt, F. C., and Hyman, J. M. (Courant Institute, New York University: preprint, 1975).
22. Smale, S., and Williams, R. (Department of Mathematics, Berkeley: preprint, 1976).
23. May, R. M., Science, 186, 645-647 (1974).
24. Moran, P. A. P., Biometrics, 6, 250-258 (1950).
25. Ricker. W. E., J. Fish. Res. Bd. Can., 11, 559-623 (1954).
26. Cook, L. M., Nature, 207, 316 (1965).
27. Macfadyen, A., Animal Ecology: Aims amid Methods (Pitman, London, 1963).
28. May, R. M., J. theor. Biol., 51, 511-524 (1975).
29. Guckenheimer, J., Proc. AMS Symposia in Pure Math., XIV, 95-124 (1970).
30. Gilbert, E. N., and Riordan, J., Illinois J. Math., 5, 657-667 (1961).
31. Preston, C. J. (King's College, Cambridge: preprint, 1976).
32. Gumowski, I., and Mira, C., C. r. hebd. Séanc. Acad. Sci., Paris, 281a, 45-48 (1975); 282a, 219-222 (1976).
33. Layzer, D., Sci. Am., 233(6), 56-69 (1975).
34. Ulam, S. M., Proc. Int. Congr. Math. 1950, Cambridge, Mass.; Vol. II, pp. 264-273 (AMS, Providence R.I., 1950).
35. Ulam, S. M., and von Neumann, J., Bull. Am. math. Soc. (abstr.), 53, 1120 (1947).
36. Kac, M., Ann. Math., 47, 33-49 (1946).
37. May, R. M., Science, 181, 1074 (1973).
38. Hassell, M. P., J. Amin. Ecol., 44, 283-296 (1974).
39. Hassell, M. P., Lawton, J. H., and May. R. M., J. Anim. Ecol. (in the press).
40. Ruelle, D., and Takens, F., Comm. math. Phys., 20, 167-192 (1971).
41. Landau, L. D., and Lifshitz, E. M., Fluid Mechanics (Pergamon, London, 1959).
42. Martin, P. C., Proc. Int. Conf. on Statistical Physics, 1975, Budapest (Hungarian Acad. Sci., Budapest, in the press).
43. Southwood, T. R. E., in Insects, Science and Society (edit. by Pimentel, D.), 151-199 (Academic, New York, 1975).
44. Metropolis, N., Stein, M. L., and Stein, P. R., Numer. Math., 10, 1-19 (1967).
45. Gumowski, I., and Mira, C., Automatica, 5, 303-317 (1969).
46. Stein, P. R., and Ulam, S. M., Rosprawy Mat., 39, 1-66 (1964).
47. Beddington, J. R., Free, C. A., and Lawton, J. H., Nature, 255, 58-60 (1975).
48. Hirsch, M. W., and Smale, S., Differential Equations, Dynamical Systems and Linear Algebra (Academic, New York, 1974).
49. Kolata, G. B., Science, 189, 984-985 (1975).
50. Smale, S. (Department of Mathematics, Berkeley: preprint, 1976).
51. May, R. M., and Leonard, W. J., SIAM J. Appl. Math., 29, 243-253 (1975).