The study of low dimensional dynamical systems which exhibit chaos is a very active area of current research. A very useful introductory account can be found in Schuster (1988). It was long thought that the complex behavior of systems of many degrees of freedom was inherently different to that of simple mechanical systems. It is now known that simple one dimensional nonlinear systems can indeed show very complex behavior. For example the family of quadratic maps fμ(x)=μx(1-x) demonstrates many of these features. This is very well described in a recent book by Devaney (1986). The connection between a discrete mapping, and the solution of a system of ordinary differential equations in a molecular dynamics simulation is clear when we realise that the numerical solution of the equations of motion for a system involves an iterative mapping of points in phase space. Although we are solving a problem which is continuous in time, the differential equation solver transforms this into a discrete time problem. The result is that if the mapping f takes Γ(0) to Γ(∆) where ∆ is the time step, then Γ(n∆)fn(Γ(0)). Here fn means the composite mapping consisting of n repeated operations of f., f(K f(Γ(0))K).
An important difference exists between difference equations and similar differential equations, for example if consider the differential equation
(10.5)
the solution can easily be obtained
(10.6)
where x0=x(t=0). The trajectory for this system is now quite straightforward to understand. The solution of the quadratic map difference equation is a much more difficult problem which is still not completely understood.
The quadratic map is defined by the equation
(10.7)
If we iterate this mapping for μ=4, starting with a random number in the interval between 0 and 1, then we obtain dramatically different behavior depending upon the initial value of x. Sometimes the values repeat; other times they do not; and usually they wander aimlessly about in the range 0 to 1. Initial values of x which are quite close together can have dramatically different iterates. This unpredictability or sensitive dependence on initial conditions is a property familiar in statistical mechanical simulations of higher dimensional systems. If we change the map to xn+1=3.839xn(1-xn) then a random initial value of x leads to a repeating cycle of three numbers (0.149888..,0.489172..,0.959299..). This mapping includes a set of initial values which behave just as unpredictably as those in the μ=4 example but due to round-off error we don't see this randomness.
Before we look at the more complicated behavior we consider some of the simpler properties of the family of quadratic maps. First we require some definitions; x1 is called a fixed point of the map f if f(x1)=x1. x1 is a periodic point, of period n, if fn(x1)=x1, where fn represents n applications of the mapping f. Clearly a fixed point is a periodic point of period one. The fixed point at x1 is stable if |f'(x1)|<1. We will consider the quadratic map fμ(x) on the interval 0<x<1, as a function of the parameter μ.
Region 1 : 0<μ<1
The mapping fμ(x) has only one fixed point x=0. f'μ(0)=μ<1 so that in this region the fixed point at x=0 is attracting (or stable).
Region 2 : 1<μ<3
fμ(x) has two fixed points x=0 and xp=(μ-1)∕μ. The fixed point x=0 is repelling (or unstable) while |f'μ(xp)|=|2-μ|<1, so that xp is an attracting (or stable) fixed point.
Region 3 :
In this region both the fixed points of fμ(x) are unstable so we consider the composite mapping f 2μ(x)=fμ(fμ(x))=μ2x(1-x)[1-μx(1-x)]. f 2μ has the fixed points of the original mapping fμ(x) at x=0 and xp, but as before both of these are unstable. f 2μ also has two new fixed points at 
. These two fixed points x± of f 2μ are points of period two in the original mapping fμ(x), (referred to as a 2-cycle). f 2μ(x±)=4+2μ-μ2 so the 2-cycle is stable for 
. If we consider finding solutions of the equation f 2μ(x)=μ2x(1-x)[1-μx(1-x)=x, then we see that we have to find the zeros of a polynomial of
order 4. This has 4 solutions; the two fixed points and x±. The 2-cycle solution x± is real for μ>3 and a complex conjugate pair for μ<3. Note however, that the two solutions x± appear at the same parameter value μ.
Region 4 ,5, etc :
The period doubling cascade where the stable 2-cycle loses its stability, and a stable 4-cycle appears; increasing μ the 4-cycle loses stability and is replaced by a stable 8-cycle; increasing μ again leads to the breakdown of the 2n-1-cycle and the emergence of a stable 2n-cycle. The μ bifurcation values get closer and closer together, and the limit as n→∞ the bifurcation value is approximately μ∞=3.5699456.
The Chaotic Region : μ∞< μ<4
Here stable periodic and chaotic regions are densely interwoven. Chaos here is characterized by sensitive dependence on the initial value x0. Close to every value of μ where there is chaos, there is a value of μ which corresponds to a stable periodic orbit, that is, the mapping also displays sensitive dependence on the parameter μ. The windows of period three, five and six are examples. From the mathematical perspective the sequence of cycles in a unimodal map is completely determined by the Sarkovskii theorem (1964). If f(x) has a point x which leads to a cycle of period p then it must have a point x' which leads to a q-cycle for every q←p where p and q are elements of the following sequence (here we read ← as precedes)
This theorem applies to values of x at a fixed parameter μ, but says nothing about the stability of the cycle or the range of parameter values for which it is observed.
Figure 10.1 The iterates of the quadratic map for some particular values of the parameter μ. The horizontal axis is xn and the vertical axis is xn+1. For μ=2 and 2.9 there is a single stable fixed point. For μ=3.3 there is a stable 2-cycle; for μ=3.5 a stable 4-cycle and for μ=3.561 a stable 8-cycle. The value μ=3.83 is in the period three window.
Figure 10.2 The iterates of the quadratic map as a function of the parameter μ. The horizontal axis is the parameter 1≤μ≤4, and the vertical axis is the iterate 0≤xn≤1.
Figure 10.3 The iterates of the quadratic map as a function of the parameter μ. This is an expanded version of Figure 10.2 to include more detail in the chaotic region. The horizontal axis is the parameter 3.5≤μ≤4, and the vertical axis is the iterate 0≤xn≤1. The windows of period three (at about μ=3.83), period five (at about μ=3.74), and period six (at about 3.63) are clearly visible.
Region ∞- : μ=4
Surprisingly for this special value of μ it is possible to solve the mapping exactly (Kadanoff, 1983).
Making the substitution 
(10.8)
A solution is 
, or 
Since xn is related to 
adding an integer to 
leads to the same value of xn. Only the fractional part of has significance. If is written in binary (base 2) notation
(10.9)
then the mapping is simply shifting the decimal point one place to the right and
removing the integer part of 
. The equivalent mapping is
(10.10)
It is easy to see that any finite precision approximation to the initial starting
value 
consisting of N digits will lose all of its significant digits in N iterations.
If x0 evolves to f(x0) after one iteration then the distribution δ(x-x0) evolves to δ(x-f(x0)) after one iteration. This can be written as
(10.11)
An arbitrary density ρn(x) constructed from a normalized sum of (perhaps infinitely many) delta functions, satisfies an equation of the form
(10.12)
The invariant measure ρ(x), or steady state distribution, is independent of time (or iteration number n) so
(10.13)
There is no unique solution to this equation as ρ(x)=δ(x-x*) where x* is an unstable fixed point of the map, is always a solution. However, in general there is a physically relevant solution and it corresponds to the one that is obtained numerically. This is because the set of unstable fixed points is measure zero in the interval [0,1] so the probability of choosing to start a numerical calculation from an unstable fixed point x*, and remaining on x*, is zero due to round off and truncation errors.
Figures 10.2 and 10.3 show the iterates of the quadratic map. In Fig. 10.4 we present the invariant measure of the quadratic map in the chaotic region. The parameter value is μ=3.65. The distribution contains a number of dominant peaks which are in fact fractional power law singularities.
Figure 10.4 The distribution function for the iterates of the quadratic map in the chaotic region, at μ=3.65. The horizontal axis is the value of the iterate, and the vertical axis is the probability. Notice the distribution of narrow peaks which dominate the probability distribution.
For the transformed mapping 
, it is easy to see that the continuous loss of information
about the initial starting point with each iteration of the map, means that the
invariant measure as a function of 
is uniform on [0.1] (that is 
). From the change of variable 
it is easy to see that x is a function of , 
(but not the reverse). If 
, then the number of counts in the distribution function
histogram bin centered at x1 with width dx1, is equal to the number of counts in the bins centered at 
and 
with widths 
That is
. (10.14)
It is then straightforward to show that the invariant measure as a function of x is given by
. (10.15)
This has inverse square root singularities at x=0 and x=1. In Figure 10.5 we present the invariant measure for the quadratic map at μ=4. The two singularities of type (x-x0)-½ at x0=0 and x0=1 are clearly shown.
Figure 10.5 The distribution of iterates for the quadratic map at μ=4. The horizontal axis is the iterate, and the vertical axis is the probability. When correctly normalized, this agrees well will equation (10.15).
Region ∞ : μ>4
Here the maximum of fμ(x) is greater than one. Once the iterate leaves the interval 0<x<1 it does not return. The mapping f 2μ(x) has two maxima, both of which are greater than one. If I is the interval [0,1], and A1 is the region of I mapped out of I by the mapping fμ(x), A2 the region of I mapped out of I by f 2μ(x), etc., then the trajectory wanders the interval defined by I-(A1∪A2∪K) It can be shown that this set is a Cantor set.
This example of a seemingly very simple iterative equation has very complex behaviour as a function of the parameter μ. As μ is increased the stable fixed point becomes unstable and is replaced by stable 2n-cycles (for n=1,2,3,K), until chaotic behaviour develops at μ∞ (about 3.5699456). For μ∞>3.5699456K the behaviour of the quadratic map shows sensitive dependence upon the parameter μ, with an infinite number of islands of periodic behaviour immersed is a sea of chaos. This system is not atypical, and a wide variety of nonlinear problems show this same behaviour. We will now consider a simple model from hydrodynamics which has had a dramatic impact in the practical limitations of weather forecasting.
Consider two flat plates, separated by a liquid layer. The lower plate is heated and the fluid is assumed to be two-dimensional and incompressible. A coupled set of nonlinear field equations must be solved in order to determine the motion of the fluid between the plates (the continuity equation, the Navier-Stokes equation and the energy equation). These equations are simplified by introducing the stream function in place of the two velocity components. Saltzman (1961) and Lorenz (1963) proceed by making the field equations dimensionless and then representing the dimensionless stream function and temperature by a spatial Fourier series (with time dependent coefficients). The resulting equations obtained by Lorenz are a three parameter family of three-dimensional ordinary differential equations which have extremely complicated numerical solutions. The equations are
(10.16)
where σ, r and b are three real positive parameters. The properties of the Lorenz equations have been reviewed by Sparrow (1982) and below we summarize the principle results.
Simple Properties
- Symmetry - The Lorenz equations are symmetric with respect to the mapping (x,y,z)→(-x,-y,z).
- The z-axis is invariant. All trajectories which start on the z-axis remain there and move toward the origin. All trajectories which rotate around the z-axis do so in a clockwise direction (when viewed from above the z=0 plane). This can be seen from the fact that if x=0, then ẋ>0 when y=0, and ẋ<0 when y<0.
- Existence of a bounded attracting set of zero volume, that is the existence of an attractor. The divergence of the flow, is given by
(10.17)
The volume element V is contracted by the flow into a volume element Vexp[-(1+b+σ)t] in time t. We can show that there is a bounded region E, such that every trajectory eventually enters E and remains there forever. There are many possible choices of Lyapunov function which describe the surface of the region E. One simple choice is V=rx2+σy2+σ(z-2r)2. Differentiating with respect to time and substituting the equations of motion gives
(10.18)
Another choice of Lyapunov function is E=r2x2+σy2+σ(z-r(r-1))2 for b≤r+1. This shows that there exists a bounded ellipsoid, and together with the negative divergence shows that there is a bounded set of zero volume within E towards which all trajectories tend.
- Fixed points. The Lorenz equations have three fixed points;
one at the origin, the other two are at
and 
. - Eigenvalues for linearized flow about the origin are
- Stability
| 0<r<1 |
The origin is stable |
| r>1 |
The origin is non-stable. Linearized flow about the origin has two negative and one positive, real eigenvalues. |
 |
C1 and C2 are stable. All three eigenvalues of the
linearized flow about C1 and C2, have negative real part. For r>1.346 (σ=10, |
 |
C1 and C2 are non-stable. Linearized flow about C1 and C2 has one negative real eigenvalue and a complex conjugate pair of eigenvalues with positive real part. |
) there is a complex conjugate pair of
eigenvalues. Again we have a nonlinear system which is well behaved for small values of the
parameter r, but for chaotic behaviour begins. Typical iterates of the Lorenz model
are shown in Fig. 10.6.
