Chaotic systems are unpredictable, but how much so? Can we say that one system is more chaotic than another? It's clear from our simulations and visualizations of chaotic attractors that they come in many shapes and forms and have distinct properties, such as being fractals and having sensitive dependence on initial conditions.
Keying off of the idea of sensitivity, here we'll explore how to measure the degree of unpredictability using Lyapunov characteristic exponents (LCEs). First, we write code for the simplest case: one-dimensional maps. Then we look at two-dimensional maps where, due to the higher dimension, there can be independent directions of instability (stretching) and stability (shrinking). (In one dimension if both occur, then they occur along one direction.) This means there are two Lyapunov exponents. From these, we'll learn how to calculate the fractal dimension of chaotic attractors and also to measure how rapidly areas contract. In turns out that all of these quantities—Lyapunov exponents, fractal dimension, and contraction rate—are related and we'll use their relationships as a check on whether our estimates of them are at least self-consistent.
The next step is to extend the Lyapunov exponents to three-dimensional flows. We'll measure the spectrum of the three Lyapunov exponents of the Lorenz chaotic attractor in three dimensions. The code, as you will see, becomes substantially more complex. The version we'll develop here is rather explicit in its calculations. We'll return in a later lab to investigate how we can use Python's object-oriented features to write more a compact and readable algorithm.
At the end, we return to some interactive visualization methods using the visual module. First, we look a 3D version of two coupled logistic maps—Rossler's folded-towel map. And, then, as a visual supplement to our measuring chaos, we use dot spreading to see how close initial conditions separate exponentially fast on the Lorenz attractor. The visual module's multi-threaded nature makes for a nice tool for dot spreading explorations.
As with the previous lab, the following notes and comments are merely a guide to the code and algorithms. The intention is that you should read and understand the programs. A number of coding techniques that will be useful in your projects are sprinkled throughout.
One-dimensional maps are the easiest case for writing a program to estimate the degree of chaos. There's a single Lyapunov exponent. And, as we saw in the theory lectures, the algorithm is fairly direct: As we iterate the map, we calculate the map's slope at each iterate, summing up the logarithm of the absolute value of the slope. Since we are interested only in the magnitude of stretching or shrinking we use the slope's absolute value.
LCE for Tent map: TentMapLCE.py: This program estimates the LCE for the tent map. The overall program structure is very similar to that for our previous 1D map simulators. The first difference is that we add a function for calculating the map's derivative at a point. The second is that we add a statement in the main map-iteration loop that accumulates the logarithm (of the absolute value) of the slope at the current state.
Recall from the theory lectures that we actually know the analytical form of how the LCE depends on the tent map control parameter. Thus, in one sense, we don't have to do the calculation numerically. But it is always a good idea when developing a new algorithm to test it out on cases for which you have analytical solutions.
Question: Does the numerical calculation verify the analytical one?
Since at each parameter value the calculations aren't too arduous we set the program up to plot the LCE as a function of the tent map's control parameter. As it loops over the parameter values, it calculates the LCE at each parameter value, storing the results in an array for later plotting (with matplotlib).
The result, though, is that this is not a fast program. We are asking it to do a fair amount of computation, especially, over a sweep of parameters. But that's what machines are supposed to do for us, right? Work that we don't want to or can't do.
There is an error waiting to occur in this program. We should have code that watches out for the logarithm being passed a zero value. Which it would certainly get if we set the tent map parameter a = 0. This, of course, should be detected. We'll come back to how to catch and respond to errors when we talk about the facilities the Python provides for this.
LCE and Bifurcation Diagram for Tent map: TentMapBifnLCE.py: Here we merge the bifurcation diagram and LCE estimation programs to produce an overall graphic that is much more informative that either alone. The bifurcation diagram does not tell us where the most stable orbits are nor that stability is lost at bifurcations nor that the degree of instability increases with parameter. In a complementary fashion, the LCE plot does not tell us what the orbit periods are nor where chaotic bands merge. Shown together we get to see these by comparing.
Note how the upper (LCE) plot x-axis tick marks are turned off and the x-axis has no label. Doing this removes clutter and so eases visual comparison between the diagrams. Though these are simple things to do, generally it's important to think carefully about the layout of a program's graphical output. Clear and simple design can make a tool much easier to use. (I recommend reading the books by Tufte in the Supplemental Reading list as an aid in designing graphical displays.)
In 2-dimensional dynamical systems we can have two directions of stretching and shrinking along an orbit. The result is that there are two LCEs and we have to monitor the separation or approach of neighboring, close points in two complementary directions.
The Jacobian matrix, which is the linearized version of the 2D map at a point, indicates how nearby points evolve forward along the orbit. We use the Jacobian, the so-called tangent space map, to evolve forward a pair of displacement vectors. This, in turn, let us monitor the divergence or convergence of nearby orbits.
The program implements the pull-back algorithm for estimating the two LCEs. The main idea is that, as we follow along an orbit, we also carry along a 2D coordinate frame attach to that orbit. The coordinate frame is represented by two vectors that we keep orthogonal and normalized. The coordinate frame vectors are evolved using the Jacobian matrix.
How much the coordinate frame vectors stretch or shrink is accumulated and become our estimates of the LCEs. Note that we must keep the vectors orthogonal so that they measure stretching in one direction and shrinking in the complementary direction. If we didn't subtract out the part of the second vector that was being stretched, then it's change in size would be dominated by stretching. What we do is subtract out the component of the first vector, which is measuring the stretching, from the second vector. The result is that the latter only points in the direction of shrinking and so its change in size is used to estimate the second LCE.
LCE for the Henon map: HenonMapLCE.py: This program implements the pull-back method for estimating the LCEs for the Henon map. Note how we implement the normalization and orthogonalization of the coordinate frame vectors. Note also that we calculate ahead of time the functional form of the Jacobian, defining a function that allows us to directly use this matrix to evolve the coordinate frame vectors.
The determinant of the Jacobian gives the area contraction rate. In the case of the Henon map the determinant is constant and so we don't have to estimate it along an orbit, which one must do in the general case.
As we discussed in the theory lectures, the area contraction rate must be the sum of the two LCEs. This makes intuitive sense. If the LCEs measure the stretching along two complementary directions, together they should measure the change in areas. We use the fact that their sum should be the area contraction rate to monitor accuracy of the LCE estimates. Generally, it's a very good idea to include these kinds of consistency checks on numerical estimates. The important thing is to find quantities that are estimated or known from independent information.
Finally, we estimate the fractal dimension from the LCEs using the Kaplan-Yorke formula described in the dynamics lectures.
Dynamical systems in three dimensions have three possible directions in which the distance to neighboring points can stretch or shrink. And so there are three LCEs. Again, we have the Jacobian matrix which monitors the local (linear) stability at a state and which we use to evolve forward displacement vectors attached to a reference orbit.
The basic approach of estimating LCEs is the same as illustrated in two dimensions. The main difference is that the dimension is larger and so there are more vectors to work with. The pull-back method generalizes in a straightforward manner, though. We have a (now 3D) coordinate frame that is evolved forward using the Jacobian, or tangent space flow.
LCE Spectrum for Lorenz ODE: LorenzODELCE.py: This program not only illustrates estimating the full spectrum of LCEs in three dimensions, but also for a continuous-time flow—the Lorenz ODEs. The main difference is that, since we're in continuous time and not discrete as above for the Henon map, we have to use the Runge-Kutta integrator for the ODEs. Additionally, this is used to evolve the state forward in time accurately along with evolving the coordinate frame vectors.
As you see the code, though a more or less straightforward generalization of the 2D discrete-time case, is much more complicated. Some of this is due to our explicitly listing out the vector calculations in component form. Mostly though, it's a more complex procedure.
Since this is in continuous time, there is one LCE that is zero, indicating a direction of neutral stability. This corresponds to the direction along orbit—a direction along which the orbit neither speeds up or slows down over the long term.
The result is that the LCE spectrum for the Lorenz attractor has one positive, one zero, and on negative LCE. That is, its spectrum has the signature of (+,0,-).
The program also measures the volume dissipation or volume contraction using the trace of the Jacobian. (We use the trace of the Jacobian matrix for this, rather than the determinant used for discrete-time maps, since the system has continuous time.) This must be the sum of the three LCEs and so we use it to monitor accuracy of LCEs estimates. We also use the closeness of the LCE associated with the direction along the orbit to zero as another measure of confidence in our estimated LCEs.
To see how bad the estimates can be, rather than running with the number of iterations and pull-backs given, try reducing these. You'll see the estimation accuracy decrease.
Let's come back to interactive 3D visualization, this time of a map in 3D.
Rossler Towel Map: TowelMap.py: Here we use mayavi to see the folded structure of Rossler's (aptly named) towel map. Look at the equations and you'll see that it's very similar to two coupled logistic maps that we studied before. This time a third coordinate has been added to reveal more of the folding structure of this chaotic map.
This program is stunningly simple; a little vacation from the complexity of the preceding programs. Of course, hidden from us is all of the interactive 3D visualization which the mayavi module provides. The simulation and visualization are very simple.
Remember that when running mayavi programs within ipython be sure to start the latter with:
ipython --pylab
The programs certainly would benefit by being rewritten to use Python's array and vector objects. As noted, this would also mean that many of the algorithms would apply to any dimension dynamical systems. They would also be closer to the underlying mathematics and so easier to understand and maintain.
They should also be broken into separate classes and modules, each of which had a clearly defined single purpose. For example, the Runge-Kutta integrator should be in its own module; as should, say, a library of different dynamical systems. There should also be a module devoted to graphical display and user interaction.
We'll return repeatedly to the issues of code complication and overall software design to manage it.