𝔉𝔲𝔯𝔪𝔞𝔫 𝔘𝔫𝔦𝔳𝔢𝔯𝔰𝔦𝔱𝔶 Electronic Journal of Undergraduate Mathematics Volume 4, 1998 1-8

QUANTIFYING CHAOS IN DYNAMICAL SYSTEMS WITH LYAPUNOV EXPONENTS

MICHAEL VAN OPSTALL

Abstract. In this paper, we analyze the dynamics of a four dimensional mechanical system which exhibits sensitive dependence on initial conditions. The aim of the paper is to introduce the basic ideas of chaos theory while assuming only a course in ordinary differential equations as a prerequisite.

1. Introduction

Dynamical systems, in short, are systems which exhibit change. As such, the field of dynamical systems is varied and rich. Many dynamical systems can be modeled by systems of differential equations or discrete difference equations. Such systems are called deterministic. Examples of such systems include those of classical mechanics. Sensitive dependence on initial conditions is a phenomenon where slight distance between the initial conditions of a system grows exponentially. Deterministic dynamical systems that exhibit a sensitive dependence on initial conditions are known as chaotic. Many physical systems are chaotic, from the driven simple pendulum to the more complex system modeled in this paper.

Dynamical systems are classified as discrete or continuous. A discrete dynamical system (given by one or more difference equations) is one in which a function f is iterated on an initial condition x_{0}. The set of all points generated by iterating f beginning with x_{0} is known as the orbit of x_{0} under f. A continuous system is generally given by one or more differential equations. Continuous orbits are known as trajectories.

There are several difficulties in working with chaotic systems. Systems of differential equations that believe chaotically are always nonlinear. This nonlinearity makes an analytic solution to these equations difficult. In addition to the nonlinearity, a continuous system which exhibits sensitive dependence on initial conditions must have dimension of at least three (that is, it must have three independent variables). The system discussed in the present paper has degree four, and hence cannot be easily visualized. Despite these difficulties, the fundamental concepts of the science are accessible to anyone who has taken a course in ordinary differential equations.

1.1. One-Dimensional Discrete Systems. Despite their simple nature, systems in a single variable can be used to model many things. One good example is the logistic map x_{n+1} = μx_{n}(l − x_{n}), which is used as a simple model for population Figure 1. The logistic map with μ = 2.9, with an attracting fixed point (left) and another with it μ = 4.0, a repelling fixed point.

Received by the editors September 2, 1998.

1991 Mathematics Subject Classification. 34D08.

Key words and phrases. Chaos, Differential Equations, Dynamical Systems.

This paper was written while the author was an undergraduate at Hope College.

growth. Some features of dynamical systems are easiest to demonstrate in singledimensional systems, so a few are described here.

The orbit of the function is computed according to the relation x_{n+1} = f(x_{n}). The logistic map described above is an example of a one-dimensional discrete system. A point where the function’s value is unaffected by further iteration (i.e. x_{n+1} = x_{n}) is called a fixed point. A fixed point which is approached by orbits is known as an attractor and one from which orbits diverge is a repeller. Figure 1 represents a logistic map with an attracting fixed point (μ = 2.9) and a chaotic logistic map with a repelling fixed point (μ = 4.0).

One way to quantify chaotic behavior in a system is to measure the divergence between orbits of two points with small initial separation. Assume f^{n} is the nth iteration of a function f. Then, for two different initial conditions, x and x + ε, the separation between these orbits is given by |f^{n}(x + ε) − f^{n}(x)|, as a function of the number of iterations. If we assume that the separation of the trajectories grows (or shrinks) exponentially we have

|f^{n}(x + ε) − f^{n}(x)| ≈ εe^{nλ},

and λ is called the Lyapunov exponent.

If we take the initial separation, ε, between trajectories to be small, we obtain

λ ≈ \frac{1}{n} log |\frac{f^{n}(x + ε) − f^{n}(x)}{ε}| ≈ \frac{1}{n} log |\frac{df^{n}}{dx} |  (1)

Noting that, x_{n} = f(x_{n−1}) = f^{n}(x), we find df^{n}/dx using the chain rule

\frac{df^{n}}{dx} = f′(f^{n−1}) · f′(f^{n−2}) ⋯ f′(x_{0}) = ∏_{m=0}^{n−1} f′(x_{m})

QUANTIFYING CHAOS IN DYNAMICAL SYSTEMS WITH LYAPUNOV EXPONENTS 3

From this, we arrive at our final formula for the Lyapunov exponent of a onedimensional discrete system:

λ = lim_{n→∞} \frac{1}{n} log | ∏_{m=0}^{n} f′(x_{m})| = lim_{n→∞} \frac{1}{n} ∑_{m=0}^{n} log |f′(x_{m})|  (2)

This exponent represents the average exponential rate of divergence of nearby orbits. A zero exponent implies linear divergence. A positive exponent indicates sensitive dependence on initial conditions, as points initially close together will diverge exponentially along neighboring trajectories. Negative exponents are found in systems where trajectories converge so the initial separation between two points will decrease in time.

Noting that the formula for a single-dimensional Lyapunov exponent is simply an average of the logarithm of the size of the derivative, a formula for continuous systems can be obtained. For a continuous system, the mean becomes the expected value of log |f′(x)| and we have the following formula for the Lyapunov exponent for a single dimensional continuous system:

λ = ∫ f(x) log |f′(x)|dx.  (3)

1.2. Higher Dimensional Systems. For many real systems, a single-dimensional model is inadequate. Unfortunately, along with a better, multi-dimensional model, we gain more problems in calculating the Lyapunov exponent of a system. The equation derived for single variable discrete systems does not directly apply, and non-linear differential equations pose problems, as they are difficult or impossible to solve. We must often resort to numerical methods to solve these problems.

1.2.1. Phase Space. The phase space of a system is the n-dimensional space in which the points of an n-dimensional system reside. A graph of trajectories in the phase space is known as a phase diagram. For two dimensional systems, the phase space lies in the plane (known as the phase plane), and is easily visualized. For higher dimensional systems, however, the phase space is often projected into two dimensions for easy viewing. In our system (described below), the four variables defining the phase space were paired to produce two phase diagrams. A twodimensional phase diagram often plots the velocity of a body against its position.

1.2.2. Attractors. In the one-dimensional case, points to which orbits converged were known as attracting fixed points. The fixed point is a special case of an attractor. In higher dimensional spaces, trajectories with small initial separation are sometimes pulled together into a single trajectory, an attractor. In these higher dimensional systems, these attractors can be curves or surfaces. An attractor in a chaotic system is known as a strange attractor.

1.3. Lyapunov Exponents. In an n-dimensional dynamical system, we have n Lyapunov exponents. Each λ_{k} represents the divergence of k-volume (k = 1: length, k = 2: area, etc.). The sign of the Lyapunov exponents indicates the behavior of nearby trajectories. A negative exponent indicates that neighboring trajectories converge to the same trajectory. A positive exponent indicates that neighboring trajectories diverge. When trajectories diverge exponentially, a slight error in measurement of the initial point could be catastrophic, as the error grows exponentially as well. If ε in equation (1) is taken to be the slight error in measuring a system’s Figure 2. Left: Convergence of trajectories (λ < 0), Center: two concentric circular trajectories (λ = 0), Right: divergence of trajectories (λ > 0).

state, eventually, this error grows in accordance with the Lyapunov exponent. Figure 2 represents the three types of trajectory behavior.

Any measurement taken has some error. The Lyapunov exponent affords us a measure of how quickly this error grows. If the Lyapunov exponent is negative, error actually decreases. Consider the damped pendulum; a slight error in measurement does not lead to a large overall error since the pendulum eventually comes to rest. We are primarily interested in systems where one (or more) of the Lyapunov exponents is positive. In accordance with our informal definition of chaos (behavior of a system exhibiting sensitive dependence on initial conditions), we can define a chaotic system as one with at least one positive Lyapunov exponent. Predictability in a system is lost here, and measurement error grows exponentially.

2. Experimental Setup

2.1. Physical Setup. The physical dynamical system we studied consisted of a pendulum of length ℓ and mass m attached to a block of mass M oscillating on the end of a spring with spring constant k. This apparatus was forced with forcing function f(t) = A cos φt + \sqrt{ℓ^{2} + A^{2} sin φt}, which is the motion of a camshaft of length ℓ displaced A units from the axis of rotation, driven with frequency φ. Sensors connected to a Realtime VAX recorded the position of the cart and angular displacement of the pendulum (x and θ, respectively).

To obtain data for the cart’s velocity and the pendulum’s angular velocity (v and ω, respectively), the data was generated by taking numerical “derivatives” (actually the slope between neighboring points). A diagram of our system appears in Figure 3.

2.2. Phase Diagrams. After recording data for different frequencies of forcing, two-dimensional phase plots were produced for v vs. x and ω vs. θ. The system was chaotic at high driving frequencies. The phase diagrams are given in Figure 4. Figure 5 illustrates chaotic and periodic time series. The problem of experimental noise is quite evident in these figures. Spurious data points can cause problems when calculating the Lyapunov exponents [2].

2.3. Equations of Motion. Equations of motion for this system can be obtained using the Lagrangian method. The Lagrangian L is as the difference, between

Figure 3. A schematic diagram of our physical system

Figure 4. Phase plots of angular velocity vs. angular position for the pendulum. The left is periodic and the right is a portion of the chaotic phase diagram.

potential and kinetic energy of a dynamical system. The position r⃗_{1}(t) of the cart and r⃗_{2}(t) of the pendulum are given by:

r⃗_{1}(t) = xi⃗
r⃗_{2}(t) = (x + ℓ sin θ)i⃗ − ℓ cos θj⃗.

The square of the velocity of each body is calculated by taking the square of the magnitude of the derivative of the position function:

∥r⃗_{1}′(t)∥^{2} = ẋ^{2}
∥r⃗_{2}′(t)∥^{2} = ẋ^{2} + 2ẋℓθ̇ cos θ + ℓ^{2}θ̇^{2}.

From the position and velocity functions, we get the Lagrangian:

L = \frac{1}{2}Mẋ^{2} + \frac{1}{2}m(ẋ^{2} + 2ẋℓθ̇ cos θ + ℓ^{2}θ̇^{2}) − \frac{1}{2}k(x − f)^{2} + mgℓ cos θ

Figure 5. Time series diagrams of angular position vs. time. The left is periodic and the right chaotic.

where f(t) is the displacement of the forcing piston at time t. The final equations of motion are derived according to the Euler-Lagrange equations \frac{d}{dt} (\frac{∂L}{∂q̇_{k}} ) − \frac{∂L}{∂q_{k}} where t he q_{k} are the coordinates of the system (x and θ in our case). These equations yield the final equations of motion for this system:

(M + m)ẍ − mℓθ̇^{2} sin θ + mℓθ̈cos θ + k(x − f) = 0  (4)
ℓθ̈ + ẍ cos θ + g sin θ = 0.  (5)

These equations can be further broken down into a system of four first-order differential equations, suitable for numerical integration.

3. Calculating Lyapunov Exponents

Given a system of differential equations, numerical integration affords us a method for determining the theoretical value for the system’s Lyapunov exponents. This method, described in detail by Wolf et al. in [2] and described below, is useful for determining positive Lyapunov exponents for chaotic systems. The system of equations of motion must be converted to strictly first order equations. For an n-dimensional system, n copies of n linearized equations are needed. This linearization is accomplished by multiplying the Jacobian matrix of partial derivatives of the n nonlinear functions by a column vector of the variables (comparable to approximating functions by a tangent line).

Each of the linearized equations determines a point in n space with a separation from the nonlinear system. We start with a sphere of states, centered on the nonlinear trajectory with linearized trajectories tangent to the sphere’s surface. This sphere is really nothing more than an orthonormal frame of vectors, but taking these vectors to form a sphere aids in visualization of the process. The initial state vectors defined in each direction are chosen to be orthonormal, each one perpendicular to the others and of unit length. Now the system is allowed to evolve over time. After a short time, the sphere of vectors has become an ellipsoid, with all vectors approaching the direction of greatest growth.

This presents a problem. If this is allowed to continue indefinitely, all the vectors will collapse onto the same vector and become indistinguishable. Additionally, if the largest vector continues to grow without limit, it will soon approach the size of the attractor (which is the metric diameter of the set of points that make up the attractor), at which point, the attractor folds back onto itself (vectors which have grown too large collapse to small vectors), causing a miscalculation (we lose the fact that the vector has grown, and perhaps note incorrectly that it has shrunk). These are the two major problems in calculating Lyapunov exponents. Both are solved simultaneously by renormalizing periodically using Gram-Schmidt orthonormalization. After a specified time, the vectors are measured, and are orthogonalized and brought back to a very small length. The process is repeated after orthonormalization so an average can be taken.

The highest Lyapunov exponent (recall that an n-dimensional system has n Lyapunov exponents) can be calculated by measuring the lengths of the largest vector n times over a period of t seconds, where ℓ_{n} is the length of the vector at measurement n. Then the Lyapunov exponent is:

λ = \frac{1}{n} ∑_{m=0}^{n−1} log \frac{ℓ_{m+1}}{ℓ_{m}} .  (6)

To find the other exponents, one must monitor the evolution of area or h-volume in the phase space. Then, using an analogous formula to the one above,replacing the lengths ℓ with area (or volume) A. The calculation then yields the sum of the first h exponents, where h is the number of dimensions of the space being measured (i.e. h = 1 represents length, h = 2 represents area.)

4. Conclusions

The calculation of Lyapunov exponents from collected data is similar to the process outlined above for differential equations, but pitfalls abound. Experimental noise is an issue, and phase space reconstruction (see [4]) must be considered for systems when less than all of the phase variables can be measured. Since trajectories are built from experimental data and not equations, no numerical integration is required, however. See [2] for greater detail on calculation of Lyapunov exponents.

Although many modern definitions of chaos are given in terms of orbits and periodicity, the definition given in the present paper is adequate for applications in many systems of classical mechanics. The aim of this paper is to encourage other projects. Examples of other chaotic systems include the driven simple pendulum, the double pendulum, or a multiple mass-spring system. These systems can be constructed in most college or university physics labs. We are grateful to Professor Paul De Young from the Hope College physics department for setting up our physical system.

In addition to observation of chaotic behavior, such a project reinforces or introduces several other skills in ordinary differential equations. Modeling using the Lagrangian or Newton’s laws can be practiced in these more complex systems. Additionally , chaos is an excellent context in which to introduce nonlinear equations and phase diagramming. Indeed, for simpler systems, it is possible to involve a computer algebra system such as Maple or Mathematica in drawing phase plots.

References

[1] G.L. Baker and J.P. Gollub, Chaotic Dynamics: An Introduction, Cambridge University Press, [2] A. Wolf, J.B. Swift., H.L. Swinney, J.A. Vastano, “Determining Lyapunov Exponents from a

Time Series,” Physica 16D (1985) 285-317.

[3] N.H. Packard, J.P. Crutchfield, J.D. Farmer, and R.S. Shaw, “Geometry from a Time Series,” Physical Review Letters (1980) 712-715.

[4] J.-C. Roux, R.H. Simoyi, H.L. Swinney, “Observation of a Strange Attractor,” Physica 8D (1983) 257-266.

[5] H.D.I. Abarbanel, R. Brown, J.J. Sidorowich, L.S. Tsimring, “The Analysis of Observed Chaotic Data in Physical Systems,” Reviews in Modern Physics (1993) 1331-1392.

[6] G. Benettin, L. Galgani, A. Giorgelli, J.-M. Strelcyn, “Lyapunov Characteristic Exponents for Smooth Dynamical Systems and for Hamiltonian Systems; A Method for Computing all of Them,” Meccanica (1980) 9-19.

[7] P. Bugl, Differential Equations: Matrices and Models, Prentice Hall, 1995.

Mathematics Department, University of Washington, Seattle, Washington, USA

   E-mail address: vnopstal@u.washington.edu

Sponsor: Timothy Pennings, Hope College, Holland, Michigan, USA

   E-mail address: pennings@math.hope.edu