Draft

STRANGE ATTRACTORS

Creating Patterns in Chaos

Julien C. Sprott

The University of Wisconsin

Madison Wisconsin

Copyright (c) 1993 by Julien C. Sprott

Contents

WHY THIS BOOK IS FOR YOU

CHAPTER 1: ORDER AND CHAOS

1.1 Predictability and Uncertainty

1.2 Bucks and Bugs

1.3 The Butterfly Effect

1.4 The Computer Artist

CHAPTER 2: WIGGLY LINES

2.1 More Knobs to Twiddle

2.2 Randomness and Pseudo-randomness

2.3 What's in a Name?

2.4 The Computer Search

2.5 Wiggles on Wiggles

2.6 Making Music

CHAPTER 3: PIECES OF PLANES

3.1 Quadratic Maps in Two Dimensions

3.2 The Butterfly Effect Revisited

3.3 Searching the Plane

3.4 The Fractal Dimension

3.5 Higher Order Disorder

3.6 Strange Attractor Planets

3.7 Designer Plaids

3.8 Strange Attractors that Don't

3.9 A New Dimension in Sound

CHAPTER 4: ATTRACTORS OF DEPTH

4.1 Projections

4.2 Shadows

4.3 Bands

4.4 Colors

4.5 Characters

4.6 Anaglyphs

4.7 Stereo Pairs

4.8 Slices

CHAPTER 5: THE FOURTH DIMENSION

5.1 Hyperspace

5.2 Projections

5.3 Other Display Techniques

5.4 Writing on the Wall

5.5 Murals and Movies

5.6 Search and Destroy

CHAPTER 6: FIELDS AND FLOWS

6.1 Beam Me up Scotty!

6.2 Professor Lorenz and Dr. Rössler

6.3 Finite Differences

6.4 Flows in Four Dimensions

6.5 Strange Attractors that Aren't

6.6 Doughnuts and Coffee Cups

CHAPTER 7: FURTHER FASCINATING FUNCTIONS

7.1 Steps and Tents

7.2 ANDs and ORs

7.3 Roots and Powers

7.4 Sines and Cosines

7.5 Webs and Wreaths

7.6 Swings and Springs

7.7 Roll Your Own

CHAPTER 8: EPILOG

8.1 How Common is Chaos?

8.2 But is it Art?

8.3 Can Computers Critique Art?

8.4 What's Left to Do?

8.5 What Good is it?

APPENDIXES

A. Annotated Bibliography

B. BASIC Program Listing

C. Other Computers and BASIC Versions

D. C Program Listing

E. Summary of Equations

F. Dictionaries of Strange Attractors

INDEX

Acknowledgments

I am indebted to Professor George Rowlands of the University of Warwick for introducing me to chaos and fractals and for countless stimulating discussions.

Professor Edward Pope of the University of Wisconsin - Madison Art Department assured me that these patterns have artistic appeal and suggested ways to display them.

Dr. Clifford Pickover of the IBM Watson Research Center the guru of computer visualization provided encouragement and suggestions during the early development of the ideas on which this book is based.

I would also like to thank Ray Valdes for carefully reading the manuscript and providing numerous helpful suggestions.

Finally I am grateful to scores of individuals who have critically viewed my attractors and who collectively raised my artistic consciousness to the point where this book could become a reality.

Why This Book Is For You

Art and science sometimes appear in juxtaposition one aesthetic the other analytical. This book bridges the two cultures. I have written it for the artist who is willing to devote a modicum of effort to understanding the mathematical world of the scientist and for the scientist who often overlooks the beauty that lurks just beneath even the simplest equations.

If you are neither artist nor scientist but own a personal computer for which you would like to find an exciting new use this book is also for you. Fractals generated by computer represent a new art form that anyone can appreciate and appropriate. You don't have to know mathematics beyond elementary algebra and you don't have to be an expert programmer. This book explains a simple new technique for generating a class of fractals called strange attractors. Unlike other books about fractals that teach you to reproduce well-known patterns this one will let you produce your own unlimited variety of displays and musical sounds with a single program. Almost none of the patterns you produce will ever have been seen before.

To get the most out of this book you will need a personal computer though it need not be a fancy one. It should have a monitor capable of displaying graphics preferably in color. Some knowledge of BASIC is useful although you can just type in the listings even if you don't understand them completely. For those of you who are C programmers I have provided an appendix with an equivalent version in C. You may find the exercises in this book an enjoyable way to hone your programming skills. As you progress through the book you will gradually develop a very sophisticated computer program. Each step is relatively simple and brings exciting new things to see and explore. Alternately you can use the accompanying disk immediately to begin making your own collection of strange attractors.

It is my hope that this book will instill in the artist a greater appreciation of science and in the scientist a greater appreciation of art and that it will bring enjoyment and satisfaction to computer enthusiasts both new and seasoned.

CHAPTER 1

Order and Chaos

This chapter lays the groundwork for everything that follows. Nearly all the essential ideas mathematical techniques and programming tools are developed here. If you master the material in this chapter the rest of the book should be smooth sailing.

1.1 Predictability and Uncertainty

The essence of science is predictability. Halley's comet will return to the vicinity of the earth in the year 2061. Not only can astronomers predict the very minute when the next solar eclipse will occur but also the best vantage point on the Earth from which to view it. Scientific theories stand or fall according to whether their predictions agree with detailed quantitative observation. Such successes are possible because most of the basic laws of nature are deterministic which means they allow us to determine exactly what will happen next from a knowledge the present conditions.

However if nature is deterministic there is no room for free will. Even human behavior would be predetermined by the arrangements of the molecules that make up our brains. Every cloud that forms or flower that grows would be a direct and inevitable result of processes set into motion eons ago and over which there is no possibility for exercising control. Perfect predictability is dull and uninteresting. Such is the philosophical dilemma that often separates the arts from the sciences.

One possible resolution was advanced in the early decades of the twentieth century when it was discovered that the quantum mechanical laws that govern the behavior of atoms and their constituents are apparently probabilistic which means they allow us to predict only the probability that something will happen. Quantum mechanics has been extremely successful in explaining the sub-microscopic world but it was never fully embraced by some including Albert Einstein who until his dying day insisted that he did not believe that God plays dice with the Universe.

Science has since the 1970's been undergoing an intellectual revolution that may be as significant as the development of quantum mechanics. It is now widely understood that deterministic is not the same as predictable. An example is the weather. The weather is governed by the atmosphere and the atmosphere obeys deterministic physical laws. However long-term weather predictions have improved very little as a result of careful detailed observations and the unleashing of vast computer resources.

The reason is that the weather exhibits extreme sensitivity to initial conditions. A tiny change in today's weather (the initial conditions) causes a larger change in tomorrow's weather and an even larger change in the next day's weather. This sensitivity to initial conditions has been dubbed the butterfly effect since it is hypothetically possible for a butterfly flapping its wings in Brazil to set off tornadoes in Texas. Since we can never know the initial conditions with perfect precision long-term prediction is impossible even when the physical laws are deterministic and exactly known. It has been shown that the predictability horizon in weather forecasting cannot be more than two or three weeks.

Unpredictable behavior of deterministic systems has been called chaos and it has captured the imagination of the scientist and non-scientist alike. The word "chaos" was introduced by Tien-Yien Li and James A. Yorke in a 1975 paper entitled "Period Three Implies Chaos." The term "strange attractors " from which this book takes its title first appeared in print in a 1971 paper entitled "On the Nature of Turbulence " by David Ruelle and Floris Takens. Some people prefer the term "chaotic attractor " since what seemed strange when first discovered in 1963 is now largely understood.

It's not hard to imagine that if a system is complicated (many springs and wheels and so forth) and hence governed by complicated mathematical equations that its behavior might be complicated and unpredictable. What has come as a surprise to most scientists is that even very simple systems described by simple equations can have chaotic solutions. However everything is not chaotic. After all we can make accurate predictions of eclipses and many other things.

An even more curious fact is that the same system can behave either predictably or chaotically depending on small changes in a single term of the equations that describe the system. For this reason chaos theory holds promise for explaining many natural processes. A stream of water for example exhibits smooth (laminar) flow when moving slowly and irregular (turbulent) flow when moving more rapidly. The transition between the two can be very abrupt. If two sticks are dropped side-by-side into a stream with laminar flow they will stay close together but if they are dropped into a turbulent stream they quickly separate.

Chaotic processes are not random; they follow rules but even simple rules can produce extreme complexity. This blend of simplicity and unpredictability also occurs in music and art. Music consisting of random notes or of an endless repetition of the same sequence of notes would be either disastrously discordant or unbearably boring. Art produced by throwing paint at a canvas from a distance or by endlessly replicating a pattern like a piece of wallpaper is similarly unlikely to have aesthetic appeal. Nature is full of visual objects such as clouds and trees and mountains as well as sounds like the cacophony of excited birds that have both structure and variety. The mathematics of chaos provides the tools for creating and describing such objects and sounds.

Chaos theory reconciles our intuitive sense of free will with the deterministic laws of nature. However it has an even deeper philosophical ramification. Not only do we have freedom to control our actions but the sensitivity to initial conditions implies that even our smallest act can drastically alter the course of history for better or for worse. Like the butterfly flapping its wings the results of our behavior are amplified with each day that passes eventually producing a completely different world than would have existed in our absence!

1.2 Bucks and Bugs

Enough philosophizing--it's time to look at a specific example. This example will require some mathematics but the equations are not difficult. The ideas and terminology are important for understanding what is to follow.

Suppose you have some money in a bank account that provides interest compounded yearly and that you don't make any deposits or withdrawals. Let's let X represent the amount of money in your account. When the time comes for the bank to credit your interest its computer does so by multiplying X by some number. If the interest rate were 10% the number would be 1.1 and your new balance would be 1.1 X. If your balance in the n'th year is Xn (where n is 1 after the first year 2 after the second and so forth) your balance in the year n +1 is

Xn +1 = R Xn (Eq. 1A)

where R is equal to 1.0 plus your interest rate. (R is 1.1 in this example.)

You probably know that such compounding leads to exponential growth. In terms of the initial amount X0 the amount in your account after n years is

Xn = X0Rn (Eq. 1B)

After 50 years at 10% yearly interest you will have $117.39 for every dollar you initially had invested. The bank can afford to do this only because of inflation and because money is loaned at an even higher interest rate.

Equation 1A is applicable to more than compound interest. It's how many of us have our salaries determined. It also describes population growth. Imagine some species of bug that lives for a season lays its eggs and then dies (thus avoiding the confusion of overlapping generations). The next year the eggs hatch and the number of bugs is some constant R times the number in the previous year. If R is less than 1 the bugs die out over a number of years and if R is greater than 1 their number grows exponentially.

You also know that exponential growth cannot go on forever whether it be bucks in the bank or people on the planet. Eventually something happens such as the depletion of resources and the growth slows and perhaps even reverses. Mass starvation disease crime and war are some of nature's mechanisms for limiting unbridled population growth. Thus we need to modify the above equation in some way if it is to model more closely growth patterns in nature.

Perhaps the simplest modification is to multiply the right-hand side of Equation 1A by a term such as (1 - X) that is equal to 1 if X is small (much less than 1) but which is less than 1 as X increases. Since the growth slows to zero and reverses as X approaches 1 we must think of X = 1 as representing some large number of dollars or bugs (say a million or a billion); otherwise we would never get very far! And so our modified equation called the logistic equation is

Xn +1 = R Xn (1 - Xn) (Eq. 1C)

Now you're going to get your first homework assignment. Take your pocket calculator and start with a small value of X say 0.1. To reduce the amount of work you have to do use a fairly large value of R say 2 corresponding to a doubling every year. Run X through Equation 1C a few times and see what happens. This process is called iteration and the successive values are called iterates. If you did it right you should see that X grows rapidly for the first couple of steps and then it levels off at a value of 0.5. The first few values should be approximately 0.1 0.18 0.2952 0.4161 0.4859 0.4996 and 0.5. Compare your results with the unbounded growth of Equation 1A.

You might have predicted the above result if you had thought to set Xn+1 equal to Xn in Equation 1C and solved for Xn. This value is called a fixed-point solution of the equation because if X ever has that value it will remain fixed there forever. It is also sometimes called a point attractor because every initial value of X between 0 and 1 is attracted to the fixed point upon repeated iteration of Equation 1C. Try initial values of X = 0.2 and X = 0.8. A fixed point is also sometimes called a critical point a singular point or a singularity.

If you're curious you might wonder what happens if you start with a value of X less than 0 such as -0.1 or greater than 1 such as 1.1. You should verify that the iterates are negative and that they get larger and larger eventually approaching minus infinity. We say that the solution is unbounded and that it attracts to infinity. Thus the values of X = 0 and X = 1 are like a watershed. Between these values the solution is bounded and outside these values it is unbounded.

The region between X = 0 and X = 1 is called a basin of attraction since it resembles a bathroom basin in which drops of water find their way to the drain from wherever they start. X = 0 is also a fixed point but it is unstable since values either slightly above or slightly below zero move away from zero. Such an unstable fixed point is sometimes called a repellor. Chaos results when two or more repellors are present; the iterates then bounce back and forth like a baseball runner caught in a squeeze play.

Equations that exhibit chaos have solutions that are unstable but bounded; the solution never settles down to a fixed value or even to a repeating pattern but neither does it move off to infinity. Sometimes we say that such equations are linearly unstable but nonlinearly stable. Small perturbations to the system grow but the growth ceases when the nonlinear terms become important as eventually they must. Another way to say it is that the fixed points are locally unstable but the system is globally stable. In such a case initial conditions are drawn to a special type of attractor called a strange attractor which is not a point or even a finite set of points but rather a complicated geometrical object whose properties constitute the subject of this book.

See what happens if you substitute X = 0 or X = 1 into the logistic equation. As a check on your calculations or in case you didn't do your homework Table 1.1 shows the successive iterates of X for each of the cases we have discussed.

Table 1.1 Iterates of the logistic equation for various initial values of X with R=2

An equation such as the logistic equation that predicts the next value of a quantity from the previous value is called an iterated map because it is like a road map in which each point on the earth is mapped to a corresponding point on a piece of paper. The logistic equation is a one-dimensional map because the various X-values can be thought of as lying along a straight line that stretches from minus infinity to plus infinity. Each iteration of the map moves every point along the line to a new position on the line. For the example above with R = 2 all the points between X = 0 and X = 1 walk toward X = 0.5 where they stop and remain. Other points run faster and faster toward the end of the line that stretches to minus infinity.

The logistic equation is an example of a quadratic iterated map so called because if you multiply out the right-hand side of Equation 1C it has not only a linear term RXn but also a quadratic (squared) term -RXn2. Quadratic maps are noninvertable because you can find Xn+1 from Xn but you can't go backwards because there are two values of Xn that produce the same Xn+1 and there is no way of knowing from which it came. For example Table 1.1 shows that X0 = 0.2 and X0 = 0.8 both produce X1 = 0.32. These are the two roots of the quadratic equation that you get if you try to solve for Xn in Equation 1C in terms of Xn+1.

The graph of Xn+1 versus Xn is a curve called a parabola. Because a parabola is not a straight line the map is said to be nonlinear. Chaos and strange attractors require a nonlinearity. The interesting and surprising behavior of nonlinear iterated maps is the basis for much of this book.

The first surprising result occurs if you iterate Equation 1C with R = 3.2 and an initial value of X in the range of 0 to 1. After a few iterations the solution will alternate between two values of approximately 0.5130 and 0.7995. This is called a period-2 limit cycle. Like the fixed point the limit cycle is another type of simple attractor. It is sometimes called a periodic or cyclic attractor.

It's not hard to see how cyclic behavior might arise in nature. If the population of beetles grows too large they will deplete the plants on whom they depend for food. With too few plants the beetles die out allowing the number of plants to recover leading to the next cycle of beetle growth and so forth.

Increase R a bit more to 3.5 and repeat the calculation. The result is a period-4 limit cycle with four values of approximately 0.5009 0.8750 0.3828 and 0.8269. If you keep increasing R by ever smaller amounts the period of the limit cycle will double repeatedly finally reaching chaotic behavior (an infinite period) at about R = 3.5699456. This value is sometimes called the Feigenbaum point after Mitchell J. Feigenbaum who discovered many of the interesting properties of one-dimensional maps.

When chaos occurs the successive iterates fluctuate in an apparently random and irreproducible manner. The chaotic behavior persists up to R = 4 except for an infinite number of small periodic windows. For R greater than 4 the solution is unbounded and the iterates attract rapidly to minus infinity.

The behavior described above can be summarized in a bifurcation diagram as shown in Figure 1-1 in which the limiting iterated values of the logistic equation after discarding the first few hundred iterates are plotted for a range of R from 2 to 4. This plot is called the Feigenbaum diagram and it resembles a tree on its side ("Feigenbaum " appropriately but coincidentally is German for "fig tree"). You see the fixed-point solution for R less than 3 the period-doubling route to chaos and the periodic windows at large R. The chaotic regions toward the right side of the figure are characterized by values of X that span a wide range and eventually fill the region densely with points.

X versus R

Figure 1-1. Bifurcation Diagram for the Logistic Equation Xn+1 = RXn (1 - Xn)

Each period doubling is called a bifurcation since a single solution splits into a pair of solutions. These splittings are called pitchfork bifurcations for obvious reasons. Note the period-3 window at about R = 3.84. The period-3 region begins abruptly when R is increased slightly from within the chaotic region to its left in what is called a tangent or saddle-node bifurcation. Careful inspection of the period-3 window shows that it also undergoes a period-doubling sequence at about R = 3.85. Solutions with every period can be found somewhere between R = 3 and R = 4.

Successive period doublings occur with ever increasing rapidity as one moves from left to right in Figure 1-1. The ratio of the width of each region to the width of the previous region approaches a constant equal to about 4.6692 called the Feigenbaum number. Even more remarkable is that this number arises in many different chaotic systems in nature as well as in the solutions of equations. The universality of the Feigenbaum number in chaos is reminiscent of the ubiquity of the number ¹ in Euclidean geometry.

With R = 4 the solutions occupy the entire interval from X = 0 to X = 1. Eventually X will take on a value arbitrarily close to any point in that interval (a characteristic called topological transitivity). Curiously however there are infinitely many initial values of X that don't lead to a chaotic solution even for R = 4. For example X0 = 0.5 and X0 = 0.75 lead to unstable fixed points while X0 = 0.345491... and X0 = 0.904508... produce an unstable period-2 limit cycle. By unstable we mean that if the initial values are wrong by even the slightest amount successive iterates will wander ever farther away.

Even through there are infinitely many non-chaotic initial values between zero and one the chance that you will find one by randomly guessing is negligible. For every such value there are infinitely many others that produce chaos. Such a seemingly paradoxical entity is an example of a Cantor set named after the nineteenth-century German mathematician Georg Cantor who is often credited with developing a mathematically rigorous concept of infinity.

A Cantor set contains infinitely many members (in fact uncountably infinitely many) but its members represent a zero fraction of the total! For example infinitely many points are required to cover completely the circumference of a circle but this number of points doesn't even begin to cover its interior. Such a collection (or set) of points although infinite in number is said to comprise a set of measure zero because the points fill a negligible portion of the plane. An attractor is a set of measure zero but its basin of attraction has a non-zero measure.

Few people would have guessed that such complexity could arise from such underlying simplicity. Furthermore the logistic equation is only the simplest of an endless variety of equations that can exhibit chaos. It is this dichotomy of simplicity and complexity that makes chaos beautiful to the mathematician and artist alike. In the bifurcation diagram of the logistic equation we have something with aesthetic appeal and it came from a simple quadratic equation!

1.3 The Butterfly Effect

If our goal is to seek chaotic behavior in the solution of equations we need a simple way to test for chaos. For this purpose we use the fact that chaotic processes exhibit extreme sensitivity to initial conditions in contrast to regular processes where different starting points usually converge to the same sequence of points on a simple attractor.

Suppose we iterate the logistic equation with two initial values of X that differ by only a tiny amount. Think of these values as representing two states of the atmosphere that differ only by the flapping of the wings of a butterfly. If successive iterates are attracted to a fixed point as they are for R = 2 the difference between the two solutions must get smaller and smaller as the fixed point is approached. A similar thing happens for a limit cycle. The difference between the two solutions will on average decrease exponentially.

If the solution is chaotic as is the logistic equation for R = 4 the successive iterates for the two cases will initially on average get farther apart; the difference usually grows exponentially. If the difference doubles on average with every iteration we say the Lyapunov exponent is 1. If it is reduced by half we say the Lyapunov exponent is -1. The name comes from the late nineteenth-century Russian mathematician Aleksandr M. Lyapunov (sometimes transliterated Liapunov or Ljapunov).

You can think of the Lyapunov exponent as the power of 2 by which the difference between two nearly equal X-values changes on average for each iteration. Thus the difference between the values changes by an average of 2L for each iteration. If L is negative the solutions approach one another; if L is positive we have sensitivity to initial conditions and hence chaos.

One way to detect chaos is to iterate the equation with two nearly equal initial values and see if after many iterations the values are closer together or farther apart. Another way is to make use of a principle of calculus that says that the difference in the solutions after one iteration divided by the difference before the iteration provided the difference is small is equal to the derivative of the equation for the map which for the logistic equation is

DXn+1 / DXn = R(1 - 2Xn) (Eq. 1D)

where DX is the difference between the two values of X. In Equation 1D ÆXn is the difference in the X values after n iterations and ÆXn+1 is the difference after n+1 iterations.

Since DX increases by the factor on the right of Equation 1D for each iteration the proper way to calculate the average is to start with a value of 1 and multiply it repeatedly by the right-hand side of Equation 1D at each iteration then divide the result by the number of iterations and finally take the logarithm to the base 2 of the absolute value of the result to get the Lyapunov exponent. If you prefer an equation the above description is equivalent to

L = S log2 |R (1 - 2Xn)| / N (Eq. 1E)

where the vertical bars mean that you are to disregard the sign of the quantity inside and S means to sum the quantity to its right from a value of n = 1 to a value of n = N where N is some large number. The larger the value of N the more accurate is the estimate of L.

Suppose you knew the value of X to within 0.01 for an iterated map with L = 1. After one iteration the uncertainty would be about 0.02 and after two iterations the uncertainty would be about 0.04 and so forth. After about seven iterations the error would exceed 1 and your prediction would be totally worthless. If the X-values are expressed as binary numbers each iteration would result in throwing away the right-most binary digit (bit). Thus the units of L are bits per iteration. The Lyapunov exponent is the rate at which information is lost when a map is iterated.

It is as if a succession of cartographers each copied maps from one another but every time one was copied it was only half as accurate as the previous one. If the original map were accurate to 1% the next copy would be accurate to 2% and the seventh generation copy would bear no relation to the original. If the Lyapunov exponent were -1 one bit of information would be gained at each iteration. Even a completely unknown initial condition would eventually be perfectly accurate as it approached the known fixed point or limit cycle. Unfortunately negative Lyapunov exponents are not the rule in cartography; otherwise all our maps would be self-correcting!

Figure 1-2 shows the Lyapunov exponent for the logistic equation using values of R from 2 to 4. The Lyapunov exponent is 1.0 at R = 4 since that value causes the interval of X from 0 to 1 to be mapped backed onto itself with a single fold at X = 0.5. Thus information is lost at a rate of 1 bit per iteration since each iterate has two possible predecessors. You can also see some of the periodic windows where L dips below zero toward the right edge of the plot. Also note that L is zero wherever a bifurcation occurs. At these points the solution is fraught with indecision over which branch to take and the initial uncertainty persists forever neither increasing nor decreasing.

L versus R

Figure 1-2. Lyapunov Exponent for the Logistic Equation

1.4 The Computer Artist

By now you have probably surmised that the operations described above are best carried out by a computer. The equations are simple but they must be applied repeatedly. This is precisely the kind of task at which computers excel.

There are dozens of computer types and programming languages to choose from. Currently the most popular computers are those based on the IBM PC running the MS-DOS or IBM-DOS operating system (hereafter simply called DOS). The most widely available programming language is BASIC (Beginner's All-purpose Symbolic Instruction Code) which usually comes bundled with the operating system software included with the computer. A version of BASIC called QBASIC has been included with DOS since version 5.0. BASIC may not be the most advanced computer language but it is one of the easiest to learn and to use its commands are close to ordinary English and it is more than adequate for our purposes. Furthermore modern versions of BASIC compare favorably with the best of the other languages.

The American National Standards Institute (ANSI) has established a standard for the BASIC language but it is somewhat limited and most versions of BASIC have many additions and embellishments. We will intentionally use a primitive dialect to ensure compatibility with most modern implementations and to simplify the translation into incompatible versions. In particular the programs in this book should run without modification under Microsoft BASICA GW-BASIC QBASIC QuickBASIC VisualBASIC for MS-DOS; Borland International TurboBASIC (no longer available); and Spectra Publishing PowerBASIC on IBM PCs or compatibles. You will be happiest using a modern compiled BASIC such as VisualBASIC or PowerBASIC on a fast computer with a math coprocessor.

Appendix C includes information on translating the computer programs into other partially incompatible dialects of BASIC as well as source code for use with VisualBASIC for Windows and Microsoft QuickBASIC for the Macintosh. Appendix D contains a translation into Microsoft QuickC. The BASIC programs use line numbers which have been obsolete since the mid 1980's but they are harmless and they provide a convenient way to reference lines of the program and to indicate where in the program a change is to be made.

If you follow sequentially through this book you will only need to add and change a few lines of the program as you meet each new idea. Your program will gradually grow more versatile as you work through the book. In the end you will have a powerful program that can reproduce all the examples in this book as well as an endless variety of new ones. Hence you should avoid the temptation to eliminate or to change the line numbers at least until you have a fully functional program. You may prefer to jump to Appendix B where you will find the complete final program which is also provided on the accompanying disk along with source listings in BASIC Microsoft QuickC Borland TurboC++ and a ready-to-run executable version of the program.

If you are an experienced programmer you might ridicule some of the quaint program listings. Many powerful programming structures such as block IF statements DO LOOPs and callable subroutines with local variables that produce beautifully structured programs are now standard but they have been avoided to allow backwards compatibility with more primitive versions of BASIC. They also often impose a small speed penalty. The dreaded GOTO statement has been used primarily to bypass blocks of code in deference to BASIC versions that don't support block IF statements. Lines of the program that are bypassed by a GOTO are usually indented. Blocks of the program contained within FOR...NEXT loops have also been indented. In the interest of structure and simplicity the programs have been written using numerous small modular subroutines each with a single entry point beginning with a comment line and a single exit point containing a RETURN statement albeit with global variables. The individual subroutines are separated with blank lines. It should be relatively easy for an experienced programmer to rewrite the program in a more modern format.

The program listing PROG01 iterates the logistic equation for R = 4 with an initial value of X = 0.05 and makes a graph of each iterate versus its predecessor. The program looks more complicated than it is because the various operations have been relegated to subroutines to provide a template for the more versatile cases to follow.

PROG01 Program for Iterating and Graphing the Logistic Equation

1000 REM LOGISTIC EQUATION

1010 DEFDBL A-Z                 'Use double precision
1030 SM% = 12                   'Assume VGA graphics
1190 GOSUB 1300                 'Initialize
1200 GOSUB 1500                 'Set parameters
1210 GOSUB 1700                 'Iterate equations
1220 GOSUB 2100                 'Display results
1230 GOSUB 2400                 'Test results
1240 ON T% GOTO 1190
1200
1210
1250 CLS
1260 END
1300 REM Initialize
1320 SCREEN SM%                 'Set graphics mode
1350 WINDOW (-.1
-.1)-(1.1
1.1)
1360 CLS
1420 RETURN
1500 REM Set parameters
1510 X = .05                    'Initial condition
1560 R = 4                      'Growth rate
1570 T% = 3
1590 LINE (-.1
-.1)-(1.1
1.1)

B
1630 RETURN
1700 REM Iterate equations
1720 XNEW = R * X * (1 - X)
2030 RETURN
2100 REM Display results
2300 PSET (X
XNEW)             'Plot point on screen
2320 RETURN
2400 REM Test results
2490 IF LEN(INKEY$) THEN T% = 0 'Respond to user key stroke
2510 X = XNEW                   'Update value of X
2550 RETURN

If when you first run the program your computer reports an error it will probably be in one of the following lines:

Line 1010: Be sure your version of BASIC supports double-precision (four-byte) floating-point variables. If it doesn't you may omit this line but you will then probably have to change the 4 in line 1560 to 3.99999 to avoid overflow resulting from round-off errors. With modern versions of BASIC and a computer with a math coprocessor there is no penalty and considerable advantage in using double precision. Because of the finite precision of computer arithmetic all cases will eventually repeat but with double precision the average number of iterations required before this happens is acceptably large.

Line 1320: Either your version of BASIC doesn't require this command or your computer or compiler doesn't support VGA graphics. Try reducing the 12 in line 1030 to a lower number until you find one that works. If none work try eliminating line 1320 altogether.

Line 1350: The WINDOW command defines the coordinates of the lower-left and upper-right corners of the graphics window for subsequent PSET and LINE commands. If your version of BASIC doesn't support this command you will have to delete this line and convert all the parameters in the PSET and LINE commands to address screen pixels. In such a case try replacing line 2300 with PSET (200 * X 200 - 200 * XNEW). One advantage of using the WINDOW command is that when a version of BASIC comes along that supports higher screen resolutions the program can be easily recompiled to take advantage of it.

Other errors: Look carefully for typographical errors or consult your BASIC manual to determine compatibility.

The correct program should produce a plot of the logistic parabola as show in Figure 1-3. Try different initial values of X (line 1510) and different values of R (line 1560) to confirm the behavior predicted for the logistic equation.

AMu%

Figure 1-3. The Logistic Parabola from PROG01

The logistic parabola comes from a chaotic solution but it doesn't look very complicated and it would hardly qualify as art. With one small change we can make things more interesting while at the same time illustrating the sensitivity to initial conditions. Instead of plotting each iterate versus its immediate predecessor we could plot it versus its second or third or fourth predecessor. Let's save the last 500 iterates and provide the option to plot X versus any one of them.

The changes that you need to make in the program PROG01 to accomplish this are shown in the listing PROG02. You can either go through the program and change or add lines as necessary or type the listing and save it in ASCII format and then use the MERGE command that is supported by many (mostly old) versions of BASIC to update the previous version of the program.

PROG02 Changes Required in PROG01 to Plot the Fifth Previous Iterate

1000 REM LOGISTIC EQUATION (5th Previous Iterate)

1020 DIM XS(499)
1040 PREV% = 5                  'Plot versus fifth previous iterate
1580 P% = 0
2210 XS(P%) = X
2220 P% = (P% + 1) MOD 500
2230 I% = (P% + 500 - PREV%) MOD 500
2300 PSET (XS(I%)
XNEW)                'Plot point on screen

With PREV% = 1 in line 1040 the result should be the same as for PROG01. However if you set PREV% equal to 2 you will see the logistic parabola change into a curve with two humps. Each time you increase PREV% by one you double the number of humps in the curve. Thus with PREV% = 5 the result as shown in Figure 1-4 has sixteen oscillations.

AMu%

Figure 1-4. The Logistic Parabola after five iterations from PROG02

Figure 1-4 provides a good graphical illustration of the sensitivity to initial conditions. The horizontal axis represents all possible initial conditions from zero to one. The vertical axis shows the value from zero to one corresponding to each initial condition after five iterations. It's not hard to see that two nearby points on the horizontal axis usually translate into two very different values along the vertical axis after five iterations. Try using PREV% = 10 and convince yourself that information about the initial condition is almost completely lost after ten iterations.

This exercise provides a good insight into the way a strange attractor is formed geometrically. The logistic parabola which began as a line (a one-dimensional object) is stretched and folded with each iteration eventually filling the entire plane (a two-dimensional object) after many iterations. Perhaps it reminds you of those taffy machines that repeatedly stretch and fold the taffy causing two nearby specks in the taffy after a while to be nowhere near one another. On average the distance between the specks initially increases at an exponential rate.

You should be able to think of many other examples of the sensitivity to initial conditions. When you stir your coffee to mix in the cream you're relying on a chaotic process. Two sticks dropped into the water close together just above a waterfall eventually end up far apart. Try laying two identical garden hoses side by side and turn on the water in each one at the same time without holding the ends. Chaotic processes are all around us. Their mathematical solutions usually produce chaotic strange attractors whose diversity and beauty we are about to explore.

CHAPTER 2

Wiggly Lines

In this chapter we will teach the computer to search for chaotic solutions of simple equations with a single variable. The solutions are segments of lines but the lines can wiggle in an incredibly complicated manner.

2.1 More Knobs to Twiddle

The logistic equation (Equation 1C) is an example of a dynamical system. Such systems are described by deterministic initial-value equations. This particular system has a single parameter R whose value determines the solution's behavior for all initial values of X within the basin of attraction. It's like a knob on a radio or on a stove that you can turn up or down to control the sound emitted by the radio or the convection in a pot of boiling soup.

There is a simple experiment you can do to observe the period-doubling route to chaos. Go into your bathroom or kitchen and turn on the tap but only slightly to produce a regular periodic pattern of drips. Now slowly open the tap until the pattern becomes chaotic. Just before the onset of chaos if you are sufficiently careful and patient you should observe one or more period doublings where the sound changes to something like--drip drip - drip drip - drip drip. The knob that controls the flow rate corresponds to the parameter R in the logistic equation. The dripping faucet has been extensively studied by Robert Shaw and discussed at length in his book The Dripping Faucet as a Model Chaotic System .

Most dynamical systems have more than one knob. Your kitchen faucet probably has independent control of the flow rate and the temperature of the water. With more knobs you might expect to increase the variety of ways the system can behave. Such knobs are called control parameters.

The formula for the most general one-dimensional quadratic iterated map is

Xn+1 = a1 + a2Xn + a3Xn2 (Eq. 2A)

Here we have three control parameters--a1 a2 and a3. By exploring all combinations of their values we expect eventually to observe every possible peculiar solution that the equation can have.

You might think that the initial condition X0 is a fourth knob but if the system is chaotic the solution is generally a strange attractor and all initial conditions within the basin of attraction look the same after many iterations. Of course there is no guarantee that a particular choice of X0 lies within the basin but values of X0 close to zero are within the basin about half the time and there are so many chaotic solutions over the range of the other three parameters that we can well afford to discard half of them.

The search for strange attractors proceeds as follows. Choose values for a1 a2 and a3 arbitrarily. Start with a value of X0 near zero. Iterate Equation 2A repeatedly until the solution either exceeds some large number in which case it is presumably unbounded or until the Lyapunov exponent becomes small or negative in which case the solution is probably a fixed point or limit cycle. In either event choose a different combination of a1 a2 and a3 and start over. If after a few thousand iterations the solution is bounded (X is not enormous) and the Lyapunov exponent is positive then it is likely that you have found a strange attractor.

2.2 Randomness and Pseudo-randomness

To choose values of a1 a2 and a3 we can use the random number generator provided with most computer languages. The random numbers thus produced are usually uniformly distributed between zero and one. You may wonder how a computer the epitome of determinism could ever produce a random number. This question deserves a digression since the answer provides yet another example of the very issues we have been discussing.

One way to produce a random number is to start with a value of X (the seed) between zero and one and iterate the logistic equation with R = 4 a hundred or so times. The result is a new number in the range of zero to one that is related to the seed in a complicated and sensitive way. This number is then used as the seed for the next random number which is produced in the same way. A given seed will produce the same sequence of random numbers but the sequence may not be the same on different computers or with different languages or even with different versions of the same language because of the way the numbers are rounded.

However this method of producing random numbers is not optimal. First the numbers are not uniformly distributed over the range. They tend to cluster near zero and one as the darkness of the right-hand side of Figure 1-1 suggests. Also multiplying a non-integer number by itself a hundred times is a relatively slow process on a computer.

Instead computers usually get their random numbers using the linear congruential method:

Xn+1 = (aXn + b) mod c (Eq. 2B)

where the mod (modulus) operation means to divide the quantity on its left (aXn + b) by the quantity on its right (c) and keep the remainder rather than the quotient. All the quantities in Equation 2B are integers. The constants a b and c are carefully chosen to maximize the number of steps required for the sequence to repeat which in any case can never exceed c. The numbers are uniformly distributed from zero to c - 1 but they can be transformed to the range zero to one by simply dividing Xn+1 by c. The numbers appear to be random but since they are produced using a deterministic procedure they are often called pseudo-random. Equation 2B is another example of a one-dimensional chaotic map called a shift map.

There are infinitely many conditions that truly random numbers should satisfy. Not only must the numbers be uniform over the interval but there should be no detectable relation between the numbers and any of their predecessors. In particular the sequence should repeat only after a very large number of steps. Most random number generators are deficient in certain ways. For example the random numbers produced by Microsoft QuickBASIC 4.5 repeat after 16 777 216 steps and this number is too small for some of our purposes.

The situation can be greatly improved by shuffling the numbers. Suppose we maintain a table of a hundred or so random numbers. When we want one we randomly take an entry from the table and replace it with a new random number. With this simple modification the pseudo-random numbers generated by the computer are sufficiently random for all our needs.

You should always remember that the sequence of random numbers generated by a digital computer will eventually repeat. You must take care to ensure that over the duration of a calculation such a repetition does not occur. You must also reseed the random number generator using a truly random seed such as one based on the time of day the program is started if you are to avoid repeating the same sequence each time you run the program.

2.3 What's in a Name?

When we begin to choose random values for the coefficients a1 a2 and a3 we are immediately confronted with two issues. The first is the range of values that the coefficients may have and the second is the amount by which two values of a coefficient must differ to produce attractors that are visibly different.

The first issue can be addressed by referring to the logistic equation (Equation 1C). When the value of R is too small (less than about 3.5) there are no chaotic solutions and when the value of R is too large (greater than 4) all the solutions are unbounded. A similar situation occurs for the more general one-dimensional quadratic map in Equation 2A. Thus we want to limit the coefficients to values whose magnitude (positive or negative) is of order unity. That is to say 0.1 is probably too small a value and 10 is probably unnecessarily large. This assumption can be verified by numerical experiment.

The second issue requires a subjective judgment of how dissimilar two attractors must look before we consider them to be different. In practice a change in one of the coefficients by an amount of order 0.1 will generally produce an object that is noticeably different. If we let each coefficient take on values ranging from -1.2 to 1.2 in steps of 0.1 we will have 25 possible values and we can associate each with a letter of the alphabet A through Y and have a convenient way to catalog and replicate the attractors. Limiting the coefficients to 25 values may seem excessively restrictive but since there are three coefficients for one-dimensional quadratic maps there are 253 or 15 625 different combinations.

The coefficients that correspond to the logistic equation with R = 4 are a1 = 0 a2 = 4 and a3 = -4 and they fall outside the range of -1.2 to 1.2. Thus for some purposes it will be convenient to take a larger range. A convenient way to extend the range is to use the ASCII (American Standard Code for Information Interchange) character set summarized in Table 2.1.

Table 2.1 ASCII Character Set and Associated Coefficient Values

Code Dec Coeff Code Dec Coeff Code Dec Coeff

32 -4.5 # 64 -1.3 ` 96 1.9

! 33 -4.4 A 65 -1.2 a 97 2.0

" 34 -4.3 B 66 -1.1 b 98 2.1

# 35 -4.2 C 67 -1.0 c 99 2.2

$ 36 -4.1 D 68 -0.9 d 100 2.3

% 37 -4.0 E 69 -0.8 e 101 2.4

& 38 -3.9 F 70 -0.7 f 102 2.5

' 39 -3.8 G 71 -0.6 g 103 2.6

( 40 -3.7 H 72 -0.5 h 104 2.7

) 41 -3.6 I 73 -0.4 i 105 2.8

* 42 -3.5 J 74 -0.3 j 106 2.9

+ 43 -3.4 K 75 -0.2 k 107 3.0

44 -3.3 L 76 -0.1 l 108 3.1

- 45 -3.2 M 77 0.0 m 109 3.2

. 46 -3.1 N 78 0.1 n 110 3.3

/ 47 -3.0 O 79 0.2 o 111 3.4

0 48 -2.9 P 80 0.3 p 112 3.5

1 49 -2.8 Q 81 0.4 q 113 3.6

2 50 -2.7 R 82 0.5 r 114 3.7

3 51 -2.6 S 83 0.6 s 115 3.8

4 52 -2.5 T 84 0.7 t 116 3.9

5 53 -2.4 U 85 0.8 u 117 4.0

6 54 -2.3 V 86 0.9 v 118 4.1

7 55 -2.2 W 87 1.0 w 119 4.2

8 56 -2.1 X 88 1.1 x 120 4.3

9 57 -2.0 Y 89 1.2 y 121 4.4

: 58 -1.9 Z 90 1.3 z 122 4.5

; 59 -1.8 [ 91 1.4 { 123 4.6

< 60 -1.7 \ 92 1.5 | 124 4.7

= 61 -1.6 ] 93 1.6 } 125 4.8

> 62 -1.5 ^ 94 1.7 ~ 126 4.9

? 63 -1.4 _ 95 1.8  127 5.0

ASCII codes from 0 to 31 are reserved for control codes--things like backspace carriage return and line feed. Codes from 128 to 255 can also be used but there is no universal character set associated with them. By making use of all the ASCII characters from 0 to 255 we can accommodate coefficients in the range of -7.7 to 17.8. The characters listed in the table will suffice for most of our needs however.

With such a coding scheme each attractor is represented by a sequence of characters with each character corresponding to one of the coefficients. The sequence can be thought of as the name of the attractor. We will preface the name with a character that indicates the type of equation we are dealing with. Let's use the letter A to represent one-dimensional quadratic maps. Thus the logistic equation coded in this way is AMu%. Note that the letters in the name are case sensitive (u and U are different) and so you should be careful when typing them. Such names may look strange which is perhaps appropriate for strange attractors and you shouldn't try to pronounce them! However they do provide a convenient and compact method for saving everything you need to reproduce an attractor.

2.4 The Computer Search

Before embarking on a search for strange attractors we need to generalize the formula for the Lyapunov exponent given in Equation 1E for the logistic equation. The generalization is easily obtained using differential calculus and the result is

L = S log2 |a2 + 2a3Xn| / N (Eq. 2C)

The program changes that are required to perform a search for strange attractors in one-dimensional quadratic iterated maps are given in the listing PROG03.

PROG03 Changes Required in PROG02 to Search for Strange Attractors in One-Dimensional Quadratic Maps

1000 REM ONE-D MAP SEARCH

1020 DIM XS(499)
A(504)
V(99)
1050 NMAX = 11000               'Maximum number of iterations
1160 RANDOMIZE TIMER            'Reseed random number generator
1360 CLS : LOCATE 13
34: PRINT "Searching..."
1560 GOSUB 2600                 'Get coefficients
1580 P% = 0: LSUM = 0: N = 0: NL = 0
1590 XMIN = 1000000!: XMAX = -XMIN
1720 XNEW = A(1) + (A(2) + A(3) * X) * X
2020 N = N + 1
2110 IF N < 100 OR N > 1000 THEN GOTO 2200
2120    IF X < XMIN THEN XMIN = X
2130    IF X > XMAX THEN XMAX = X
2140    YMIN = XMIN: YMAX = XMAX
2200 IF N = 1000 THEN GOSUB 3100        'Resize the screen
2250 IF N < 1000 OR XS(I%) <= XL OR XS(I%) >= XH OR XNEW <= XL OR XNEW >= XH THEN GOTO 2320
2410 IF ABS(XNEW) > 1000000! THEN T% = 2    'Unbounded
2430 GOSUB 2900                 'Calculate Lyapunov exponent
2460 IF N >= NMAX THEN T% = 2   'Strange attractor found
2470 IF ABS(XNEW - X) < .000001 THEN T% = 2 'Fixed point
2480 IF N > 100 AND L < .005 THEN T% = 2    'Limit cycle
2600 REM Get coefficients
2660 CODE$ = "A"
2680 M% = 3
2690    FOR I% = 1 TO M%        'Construct CODE$
2700        GOSUB 2800          'Shuffle random numbers
2710        CODE$ = CODE$ + CHR$(65 + INT(25 * RAN))
2720    NEXT I%
2730 FOR I% = 1 TO M%           'Convert CODE$ to coefficient values
2740    A(I%) = (ASC(MID$(CODE$
I% + 1
1)) - 77) / 10
2750 NEXT I%
2760 RETURN
2800 REM Shuffle random numbers
2810 IF V(0) = 0 THEN FOR J% = 0 TO 99: V(J%) = RND: NEXT J%
2820 J% = INT(100 * RAN)
2830 RAN = V(J%)
2840 V(J%) = RND
2850 RETURN
2900 REM Calculate Lyapunov exponent
2910 DF = ABS(A(2) + 2 * A(3) * X)
3030 IF DF > 0 THEN LSUM = LSUM + LOG(DF): NL = NL + 1
3040 L = .721347 * LSUM / NL
3070 RETURN
3100 REM Resize the screen
3120 IF XMAX - XMIN < .000001 THEN XMIN = XMIN - .0000005: XMAX = XMAX + .0000005
3130 IF YMAX - YMIN < .000001 THEN YMIN = YMIN - .0000005: YMAX = YMAX + .0000005
3160 MX = .1 * (XMAX - XMIN): MY = .1 * (YMAX - YMIN)
3170 XL = XMIN - MX: XH = XMAX + MX: YL = YMIN - MY: YH = YMAX + MY
3180 WINDOW (XL
YL)-(XH
YH): CLS
3310 LINE (XL
YL)-(XH
YH)

B
3460 RETURN

There are several things to note about the listing PROG03:

1. The maximum number of iterations (NMAX in line 1050) has been set arbitrarily to 11 000. This is the number of iterations after which a strange attractor is assumed to have been found if the magnitude of X is less than one million and the Lyapunov exponent is positive (actually greater than 0.005). You may wish to decrease NMAX to speed the rate at which attractors are found or increase NMAX if you have a very fast computer or want to give the displays more time to develop. The number of iterations is a parameter that you can adjust for the most visually appealing result. Most of the figures in this book were made with NMAX between about half a million and ten million and required between about a minute and an hour to produce.

2. The seed for the random number generator is taken in line 1160 as the number of seconds lapsed since midnight (TIMER). This choice ensures that a new sequence of random numbers is produced each time the program is run except in the unlikely event that it is run at exactly the same time each day.

3. After 1000 iterations (line 2200) the screen is resized and erased by the subroutine in lines 3100 through 3460 using the minimum and maximum values of X between the hundredth and thousandth iteration allowing a ten percent border around the attractor.

4. To save time the difference between each value of X and its predecessor is evaluated in line 2470 and if the difference is less than one millionth the solution is assumed to be a fixed point even if the Lyapunov exponent is still positive.

5. The Lyapunov exponent is not used as a criterion until after 100 iterations (line 2480) to ensure that its value is reasonably accurate.

6. The coefficients of the equation are chosen in line 2710 using random numbers which have been shuffled by the subroutine in lines 2800 though 2850 to minimize the chance of repeating the same search sequence.

The criterion for detecting a strange attractor is somewhat subjective. There will always be borderline cases for which no amount of computing will suffice to distinguish between a strange attractor and a periodic solution with a very long period. However our interest here is in finding visually interesting attractors quickly and so we can afford to make occasional mistakes. Such mistakes account for only a small number of cases.

Of the 15 625 combinations of coefficients exactly 364 (2.3%) of them are chaotic by these criteria. Some of the more visually interesting ones are shown in Figures 2-1 through 2-4 in which the values are plotted versus their fifth previous iterate. For each case the code and the Lyapunov exponent are shown at the top of the graph.

AXBH

Figure 2-1. One-Dimensional Quadratic Map

ABDU

Figure 2-2. One-Dimensional Quadratic Map

ACAV

Figure 2-3. One-Dimensional Quadratic Map

AXDA

Figure 2-4. One-Dimensional Quadratic Map

The search for strange attractors is potentially time-consuming if you have an old computer without a math coprocessor or if you are using a BASIC interpreter rather than a compiler. Even if the search is reasonably fast on your computer be forewarned that it will slow down considerably as you advance to the more complicated equations later in the book. Perhaps this is a good time to summarize some of your options for making the program run faster.

When comparing calculation speeds of various computers and compilers it is important to do the comparison with the actual program or with a benchmark that accurately reflects its mix of instructions graphics and disk access. With computer speeds doubling approximately every two years speed will eventually cease to be a consideration for the calculations described in this book. Meanwhile you need to consider the alternatives.

Table 2.2 lists the average number of strange attractors found by PROG03 per hour using various versions of BASIC on a 33 MHz 80486DX-based computer with and without a math coprocessor. The exact numbers are less important than the relative values. They provide a good indication of how the various versions of BASIC compare on calculations of the type that will be used throughout this book.

Table 2.2 Strange Attractors Found per Hour by PROG03 with Various Versions of BASIC

QuickBASIC and VisualBASIC for MS-DOS can be run from the editor environment where they function much like an interpreter or they can be used to compile a stand-alone executable program. VisualBASIC can be compiled with either of two floating point math packages; the alternate package is faster for machines without a coprocessor and the emulate package is faster for machines with a coprocessor. TurboBASIC is now obsolete and has been replaced by PowerBASIC. PowerBASIC like VisualBASIC can be compiled with either of two floating point math packages; the procedure package is similar to the VisualBASIC alternate package. A third math package NPX (87) is the same as emulate except that it will not work on a machine without a math coprocessor. The tests were done with all error trapping turned off which is inadvisable until you have a thoroughly debugged program.

If you launch the program from Microsoft Windows you might find the computation speeds considerably different from those in Table 2.2. In one test the PowerBASIC speeds were cut in half and the QuickBASIC speeds were increased slightly from the values obtained when the program was run directly from DOS. You should do your own speed tests to see what configuration provides the optimum performance on your computer and operating system.

The executable program on the disk that accompanies this book was compiled with PowerBASIC using the procedure package. If you have PowerBASIC and a math coprocessor you can recompile the program using the emulate or NPX (87) package to achieve a slight improvement in speed.

2.5 Wiggles on Wiggles

The preceding figures consist of segments of wiggly lines and thus they are not very artistic. To make things more interesting we can consider one-dimensional maps of higher order. By this we mean that we will not stop with quadratic (X2) maps but we will consider equations containing cubic (X3) quartic (X4) quintic (X5) and even higher terms.

In one sense considering higher-order terms is equivalent to plotting each iterate versus an iterate earlier than the immediately previous one. For example two successive iterations of the second-order Equation 2A yields

Xn+2 = a1(1+a2+a1a3) + (a3a2+2a1a3)Xn

(Eq. 2D)

+ a3(a2+2a1a3+a22)Xn2 + 2a2a32Xn3 + a33Xn4

which is a fourth-order polynomial. However there are only three parameters a1 a2 and a3 from which the five coefficients are uniquely determined.

A simpler and more general procedure is to allow each term in the polynomial to have its own coefficient which for fifth order gives

Xn+1 = a1 + a2Xn + a3Xn2 + a4Xn3 + a5Xn4 + a6Xn5 (Eq. 2E)

With six coefficients each with 25 possible values there are 256 or about 244 million different combinations. Even if only a small percentage of them is chaotic we would have to look at one every second for about a year before we would see them all.

The generalization of the expression for the Lyapunov exponent for a fifth-order map is given by

L = S log2 |a2 + 2a3Xn + 3a4Xn2 + 4a5Xn3 + 5a6Xn4| / N (Eq. 2C)

With these equations in hand you can easily modify the program as in PROG04 to search for one-dimensional attractors of up to fifth order. In our coding scheme a first letter of B will represent third order C will represent fourth order and D will represent fifth order. The program is written so that even higher orders can be produced by changing the quantity OMAX% in line 1060.

PROG04 Changes Required in PROG03 to Search for Strange Attractors in One-Dimensional Maps of order up to OMAX%

1000 REM ONE-D MAP SEARCH (Polynomials up to 5th Order)

1060 OMAX% = 5                  'Maximum order of polynomial
1720 XNEW = A(O% + 1)
1730 FOR I% = O% TO 1 STEP -1
1830    XNEW = A(I%) + XNEW * X
1930 NEXT I%
2650 O% = 2 + INT((OMAX% - 1) * RND)
2660 CODE$ = CHR$(63 + O%)
2680 M% = O% + 1
2910 DF = 0
2930 FOR I% = O% TO 1 STEP -1
2940    DF = I% * A(I% + 1) + DF * X
2970 NEXT I%
3000 DF = ABS(DF)

The listing in PROG04 produces a more interesting array of shapes samples of which are shown in Figures 2-5 through 2-10. The objects are still segments of lines but the wiggles themselves have wiggles and the underlying determinism is less obvious than before.

BZEZK

Figure 2-5. One-Dimensional Cubic Map

CBLCTX

Figure 2-6. One-Dimensional Quartic Map

CUTXJE

Figure 2-7. One-Dimensional Quartic Map

DBOGIZI

Figure 2-8. One-Dimensional Quintic Map

DFBIEVV

Figure 2-9. One-Dimensional Quintic Map

DOOYRIL

Figure 2-10. One-Dimensional Quintic Map

2.6 Making Music

If the preceding figures don't qualify as art perhaps they qualify as music. Since the quantity X behaves in a deterministic yet unpredictable way it may be that a sequence of musical notes determined by X will mimic the order and unpredictability that characterize music. It's easy to test.

Suppose we allow the notes to span three octaves from A-220 to A-1760. The letter refers to the musical note and the numbers refer to the frequency in cycles per second (called Hertz). We'll allow the notes to take one of twelve distinct values corresponding to the even-tempered scale and for simplicity take all the notes to be of the same duration. Thus the range of possible values of X is divided into 36 intervals and each successive iterate of X is converted into the corresponding musical note. The listing PROG05 shows the changes necessary to accomplish this.

PROG05 Changes Required in PROG04 to Produce Chaotic Music

1000 REM ONE-D MAP SEARCH (With Sound)

1100 SND% = 1                   'Turn sound on
2310    IF SND% = 1 THEN GOSUB 3500     'Produce sound
2490 Q$ = INKEY$: IF LEN(Q$) THEN GOSUB 3600        'Respond to user command
3500 REM Produce sound
3510 FREQ% = 220 * 2 ^ (CINT(36 * (XNEW - XL) / (XH - XL)) / 12)
3520 DUR = 1
3540 SOUND FREQ%
DUR: IF PLAY(0) THEN PLAY "MF"
3550 RETURN
3600 REM Respond to user command
3610 T% = 0
3630 IF ASC(Q$) > 96 THEN Q$ = CHR$(ASC(Q$) - 32)
3770 IF Q$ = "S" THEN SND% = (SND% + 1) MOD 2: T% = 3
3800 RETURN

The program allows you to toggle the sound on and off by pressing the S key. Pressing any other key exits the program. You might wish to experiment with the duration DUR of the SOUND statement in line 3520. Increasing its value from 1 (corresponding to approximately 0.055 seconds) makes the sounds more musical but the calculation will then take longer.

The use of sound to help interpret data generated by a computer is a technique that is relatively unexplored. The method is sometimes called sonification. In some cases patterns and structure in data can be more readily discerned audibly than visually. This technique was used to advantage in interpreting data from the Voyager spacecraft as it detected plasma waves near Jupiter and micrometeorites as it crossed through the rings of Saturn. The repetitive sound of a simple limit cycle contrasts sharply with the non-repetitive waverings of a chaotic time series.

CHAPTER 3

Pieces of Planes

Whereas the last chapter dealt with one-dimensional maps whose graphs are segments of lines this chapter will deal with two-dimensional maps whose graphs are pieces of planes and which thus produce much more interesting displays. This chapter provides the minimum tools for creating attractors that genuinely qualify as art. With only the information contained here the variety of available patterns is so great that you hardly need to proceed beyond this chapter but if you do stop here you will miss some delightful surprises.

3.1 Quadratic Maps in Two Dimensions

In the discussion so far the maps have involved a single variable X whose value changes with each iteration of the equation. Such maps are said to be one-dimensional because the values of X can be thought of as lying along a line and a line is a one-dimensional object. By plotting each value of X versus a previous value of X the line can be made to wiggle with considerable complexity but it always remains a line and lines are of limited interest and beauty.

The situation is more interesting if we consider iterated maps that involve two variables X and Y. In such a case each iterate produces a point in a plane where X by convention represents the horizontal coordinate and Y represents the vertical coordinate of the point. With successive iteration the points fill in some portion of the plane. The visually interesting cases as usual are the chaotic ones.

Such two-dimensional maps might arise for example from an ecological model only slightly more complicated than the logistic equation. A classic example is the predator-prey problem in which X represents the prey and Y the predator. In a simple linear model the solution is a fixed point (a unique number of rabbits and foxes) or a limit cycle (the number of rabbits and foxes oscillate reaching their maximum values at different times but eventually repeating). When nonlinear terms are introduced into the model the population of each species can behave chaotically. You can think of each point that makes up such an attractor as the population of rabbits and foxes in successive years. Since such complexity arises from these very simple models it's easy to understand why ecologists might have trouble predicting the fate of biological species!

Perhaps the best known chaotic two-dimensional map is the Hénon map (proposed by the French astronomer Michel Hénon in 1976) whose equations are

Xn+1 = 1 + aXn2 + bYn

(Eq. 3A)

Yn+1 = Xn

The quantities a and b are the control parameters analogous to R in the logistic equation. Hénon used the values a = -1.4 and b = 0.3. The necessary nonlinearity is provided by the X2-term in the first equation. The Hénon map is special in that the net contraction of a set of initial points covering an area of the XY-plane is constant with each iteration. The area occupied by the points is 30% of the area at the previous iteration (from the bYn-term). Other values of b can be used but not all values produce chaotic solutions. Unlike the logistic map the Hénon map is invertable; there is a unique value for Xn and Yn corresponding to each Xn+1 and Yn+1. You may have seen an alternate form of the Hénon equations in which the factor b appears instead in the second equation. The result of repeated iteration of Equations 3A is shown in Figure 3-1.

EWM?MPMM

Figure 3-1. The Hénon Map

The resulting graph is more than a line but less than a surface. What resembles a single line is a pair of lines each one of which in turn is another pair of lines and so forth to however close you look or whatever magnification you choose. This self-similarity is a common characteristic of a class of objects that are called fractals.

Fractals are to chaos what geometry is to algebra--the visual expression of the mathematical idea. Approaching an understanding of chaos through such visual means is appealing to those with an aversion to conventional mathematics. The Euclidean geometry that we learned in high school originated with the ancient Greeks and was developed more fully by the French mathematician Descartes and others in the 1600's. It deals with simple shapes such as lines circles and spheres. Euclidean geometry is now being augmented by fractal geometry whose father and champion is the contemporary mathematician Benoit Mandelbrot. Fractals appeared in art such as in the drawings of the Dutch artist Maurits C. Escher before they were widely appreciated by mathematicians and scientists.

Some fractals are exactly self-similar which means that they look the same no matter how much you magnify them. Others such as most of the ones in this book only have regions that are self-similar. There is no part of the Hénon attractor where you can zoom in and find a miniature replica of the entire attractor. Other fractals are only statistically self-similar which means that a magnified portion of the object has the same amount of detail as the whole but it is not an exact replica of it. Nearly all strange attractors are fractals but not all fractals arise from strange attractors.

The Hénon map produces an object with a fractal dimension that is a fraction intermediate between one and two. The fractal dimension is a useful quantity for characterizing strange attractors. Isolated points have dimension zero line segments have dimension one surfaces have dimension two and solids have dimension three. Strange attractors generally have non-integer dimensions.

Since the Hénon map has X2 as its highest-order term it is a quadratic map. The most general two-dimensional quadratic iterated map is

Xn+1 = a1 + a2Xn + a3Xn2 + a4XnYn + a5Yn + a6Yn2

(Eq. 3B)

Yn+1 = a7 + a8Xn + a9Xn2 + a10XnYn + a11Yn + a12Yn2

Equations 3B have twelve coefficients. For the Hénon map a1 = 1 a3 = -1.4 a5 = 0.3 a8 = 1 and the other coefficients are zero. If we use the initial letter E to represent two-dimensional quadratic maps the code for the Hénon map according to Table 2.1 is EWM?MPM2WM4 where we have introduced the shorthand M2 for MM and M4 for MMMM.

Values of a in the range of -1.2 to 1.2 are sufficient to produce an enormous variety of strange attractors. With increments of 0.1 there are 2512 or about 6 x 1016 different cases of which approximately 1.6% or about 1015 are chaotic. Viewing them all at a rate of one per second would require over 30 million years! Stated differently if each one were printed on an 8 1/2 by 11-inch sheet of paper the collection would cover nearly the entire land mass of the Earth.

Note that not all the cases are strictly distinct. For example if you replace X with Y and Y with -X in Equations 3B you will produce an attractor rotated 90 degrees counter-clockwise from the original. When you do this be sure to change Xn+1 and Yn+1 as well as Xn and Yn. Thus the code EM4CMWJM3? will produce a rotated version of the Hénon map. In the same fashion you can rotate an attractor through 180 degrees by replacing X with -X and Y with -Y and through 270 degrees by replacing X with -Y and Y with X . Perhaps it's easier just to rotate your computer monitor!

Besides rotations there are cases that correspond to reflections. When viewed in a mirror the attractors have left and right reversed but up and down remain the same. A transformation in which X is replaced with -X will accomplish this. Thus the code for a reflected Hénon map is ECM[MJM2CM4. In addition the reflections can be rotated. Thus there are at least eight so-called degenerate states for each attractor corresponding to rotations and reflections. Such symmetries and degeneracies play an important role in science; they often reduce the amount of work we have to do and provide relations between phenomena that initially appear different.

There are additional degenerate cases corresponding to scale changes. For example if you replace X by mX and Y by nY with m = n the attractor remains the same except it is reduced in size by a factor of m. Some of the coefficients are likely to be outside the allowed range however. The Hénon map with m = n = 2 can be generated with the code ERM1MPM2WM4. With m not equal to n the horizontal and vertical dimensions are scaled differently but since the computer rescales the attractor to fit the screen the visual result is the same.

These degeneracies show that there are many ways to code a particular attractor. Although this is true there are so many different possible combinations of coefficients that it is very unlikely that two degenerate cases will be found spontaneously. Thus the examples displayed in this chapter represent but a tiny fraction of the possibilities and you will be generating many other cases almost none of which have been seen before.

3.2 The Butterfly Effect Revisited

Two-dimensional chaotic iterated maps also exhibit sensitivity to initial conditions but the situation is more complicated than for one-dimensional maps. Imagine a collection of initial conditions filling a small circular region of the XY-plane. After one iteration the points will have moved to a new position in the plane but they will now occupy an elongated region called an ellipse. The circle will have contracted in one direction and expanded in the other. With each iteration the ellipse will get longer and narrower eventually stretching out into a long filament. The orientation of the filament also changes with each iteration and it wraps up like a ball of taffy.

Thus two-dimensional chaotic maps have not a single Lyapunov exponent but two--a positive one corresponding to the direction of expansion and a negative one corresponding to the direction of contraction. The signature of chaos is that at least one of the Lyapunov exponents is positive. Furthermore the magnitude of the negative exponent has to be greater than the positive one so that initial conditions scattered throughout the basin of attraction contract onto an attractor that occupies a negligible portion of the plane. The area of the ellipse continually decreases even as it stretches to an infinite length.

There is a proper way to calculate both of the Lyapunov exponents. For the mathematically inclined the procedure involves summing the logarithms of the eigenvalues of the Jacobian matrix of the linearized transformation with occasional Gram-Schmidt reorthonormalization. This method is slightly complicated and so we will instead devise a simpler procedure sufficient for determining the largest Lyapunov exponent which is all we need to test for chaos.

Suppose we take two arbitrary but nearby initial conditions. The first few iterations of the map may cause the points to get closer together or farther apart depending upon the initial orientation of the two points. Eventually the points will come arbitrarily close in the direction of the contraction but they will continue to separate in the direction of the expansion. Thus if we wait long enough the rate of separation will be governed only by the largest Lyapunov exponent. Fortunately this usually takes just a few iterations.

However because the separation grows exponentially for a chaotic system the points quickly become too far apart for an accurate estimate of the exponent. This problem can be remedied if after each iteration the points are moved back to their original separation along the direction of the separation. The Lyapunov exponent is then determined by the average of the distance they must be moved for each iteration to maintain a constant small separation. If the two solutions are separated by a distance dn after the n'th iteration and the separation after the next iteration is dn+1 the Lyapunov exponent is determined from

L = S log2 (dn+1 / dn) / N (Eq. 3C)

where the sum is taken over all iterations from n = 0 to n = N-1. After each iteration the value of one of the iterates is changed to make dn+1 = dn. For the cases here dn is taken equal to 10-6. This procedure will also allow us to deal with maps of three and even higher dimension in which there are more than two Lyapunov exponents.

3.3 Searching the Plane

We now have all the tools in hand to conduct a computer search for attractors in two dimensions. The procedure is the same as for one-dimensional maps except that the Lyapunov exponent calculation is done differently and the X and Y variables are plotted as a point in the plane after each iteration. Listing PROG06 shows the changes needed to accomplish such a search.

PROG06 Changes Required in PROG05 to Search for Two-Dimensional Quadratic Strange Attractors

1000 REM TWO-D MAP SEARCH

1060 OMAX% = 2                  'Maximum order of polynomial
1070 D% = 2                     'Dimension of system
1100 SND% = 0                   'Turn sound off
1520 Y = .05
1550 XE = X + .000001: YE = Y
1590 XMIN = 1000000!: XMAX = -XMIN: YMIN = XMIN: YMAX = XMAX
1720 XNEW = A(1) + X * (A(2) + A(3) * X + A(4) * Y)
1730 XNEW = XNEW + Y * (A(5) + A(6) * Y)
1830 YNEW = A(7) + X * (A(8) + A(9) * X + A(10) * Y)
1930 YNEW = YNEW + Y * (A(11) + A(12) * Y)
2140    IF Y < YMIN THEN YMIN = Y
2150    IF Y > YMAX THEN YMAX = Y
2240 IF D% = 1 THEN XP = XS(I%): YP = XNEW ELSE XP = X: YP = Y
2250 IF N < 1000 OR XP <= XL OR XP >= XH OR YP <= YL OR YP >= YH THEN GOTO 2320
2300    PSET (XP
YP)                   'Plot point on screen
2410 IF ABS(XNEW) + ABS(YNEW) > 1000000! THEN T% = 2    'Unbounded
2470 IF ABS(XNEW - X) + ABS(YNEW - Y) < .000001 THEN T% = 2
2520 Y = YNEW
2660 CODE$ = CHR$(59 + 4 * D% + O%)
2680 M% = 1: FOR I% = 1 TO D%: M% = M% * (O% + I%): NEXT I%
2910 XSAVE = XNEW: YSAVE = YNEW: X = XE: Y = YE: N = N - 1
2930 GOSUB 1700                 'Reiterate equations
2940 DLX = XNEW - XSAVE: DLY = YNEW - YSAVE
2960 DL2 = DLX * DLX + DLY * DLY
2970 IF CSNG(DL2) <= 0 THEN GOTO 3070       'Don't divide by zero
2980    DF = 1000000000000# * DL2
2990    RS = 1 / SQR(DF)
3000    XE = XSAVE + RS * (XNEW - XSAVE): YE = YSAVE + RS * (YNEW - YSAVE)
3020    XNEW = XSAVE: YNEW = YSAVE
3030    IF DF > 0 THEN LSUM = LSUM + LOG(DF): NL = NL + 1
3040    L = .721347 * LSUM / NL
3110 IF D% = 1 THEN YMIN = XMIN: YMAX = XMAX

This program produces an incredible variety of interesting patterns a small selection of which is shown in Figures 3-2 through 3-17. You will want to admire the beauty and variety of these figures and then make some of you own by running the program and if you computer has a printer using the Print Screen key to print any that you find especially appealing.

EAGHNFOD

Figure 3-2. Two-Dimensional Quadratic Map

EBCQAFMF

Figure 3-3. Two-Dimensional Quadratic Map

EDSYUECI

Figure 3-4. Two-Dimensional Quadratic Map

EELXAPXM

Figure 3-5. Two-Dimensional Quadratic Map

EEYYMKTU

Figure 3-6. Two-Dimensional Quadratic Map

EJTTSMBO

Figure 3-7. Two-Dimensional Quadratic Map

ENNMJRCT

Figure 3-8. Two-Dimensional Quadratic Map

EOUGFJKD

Figure 3-9. Two-Dimensional Quadratic Map

EQKOCSID

Figure 3-10. Two-Dimensional Quadratic Map

EQLOIARX

Figure 3-11. Two-Dimensional Quadratic Map

ETJUBWED

Figure 3-12. Two-Dimensional Quadratic Map

ETSILUND

Figure 3-13. Two-Dimensional Quadratic Map

EUEBJLCD

Figure 3-14. Two-Dimensional Quadratic Map

EVDUOTLR

Figure 3-15. Two-Dimensional Quadratic Map

EWLKWPSM

Figure 3-16. Two-Dimensional Quadratic Map

EZPMSGCN

Figure 3-17. Two-Dimensional Quadratic Map

If you are an experienced programmer you might consider writing a screen-saver program based on PROG06. Such a terminate-and-stay-resident (TSR) program is run once when the computer is turned on and leaves a portion of itself in memory constantly monitoring keyboard and mouse activity. When there is no user activity for five minutes say it would blank the screen and begin displaying a succession of unique strange attractors to prevent screen burn-in. The original screen would be restored whenever a key is pressed or the mouse is moved. PowerBASIC version 3.0 allows you to do this easily by embedding the program within POPUP statements.

3.4 The Fractal Dimension

The previous figures differ considerably in how densely they fill the plane. Some are very thin; others are thick. A good contrast is provided by Figures 3-16 and 3-17. Figure 3-16 resembles a piece of string that has been laid down in a complicated shape on the page whereas Figure 3-17 looks like a twisted piece of paper with many holes in it. Thus the object in Figure 3-16 has a fractal dimension close to 1 and the object in Figure 3-17 has a fractal dimension closer to 2.

It is possible to be more explicit and to calculate the fractal dimension exactly. Consider two simple cases one in which successive iterates lie uniformly along a straight line that goes diagonally across the page and the other in which successive iterates gradually fill the entire plane as if they were grains of pepper sprinkled on the paper from a great height. The first case has dimension 1 and the second has dimension 2. How would we calculate the dimension given the X and Y coordinates of an arbitrary collection of such points?

One method would be to draw a small circle somewhere on the plane that surrounds at least one of the points. Then draw a second circle with the same center but with twice the radius. Now count the number of points inside each circle. Let's say the smaller circle encloses N1 points and the larger circle encloses N2 points. Obviously N2 is greater than or equal to N1 since all the points inside the inner circle are also inside the outer circle.

If the points are widely separated then N2 will equal N1. If the points are part of a straight line the larger circle will on average enclose twice as many points as the smaller circle but if the points are part of a plane the larger circle will on average enclose four times as many points as the smaller circle since the area of a circle is proportional to the square of its radius. Thus for these simple cases the dimension is given by

F = log2 (N2 / N1) (Eq. 3D)

It is hardly surprising that if you do this operation for the cases shown in the figures the quantity F is neither 0 nor 1 nor 2 but rather a fraction.

With real data a number of practical considerations determine the accuracy of the result and the amount of computation required to obtain it:

1. Is a circle the best shape or would a square rectangle triangle or some other shape be better?

2. How large should the circle be?

3. Is doubling the size of the circle optimal or would some other factor be better?

4. Where should the circles be placed and how many circles are required to obtain a representative average?

5. How many points are needed to produce a reliable fractal dimension?

Let's address each of these questions in turn.

There is nothing special about circles. A rectangle triangle or any other two-dimensional figure would suffice since the area scales as the square of the linear dimension in each case. However a circle is convenient since it is easy to tell whether a given point is in its interior by comparing the radius of the circle with the distance of the point from its center.

The optimum size of the circle represents a compromise. Ideally the circles should be invisibly small since the dimension is defined in the limit of infinite resolution. However if the circles are too small they contain too few points to produce a statistically meaningful result unless an unreasonably large number of iterations is performed. We will somewhat arbitrarily take the circles to have a radius equal to about 2% of the diagonal of the smallest rectangle that contains the attractor.

Similarly doubling the radius of the circle is arbitrary. Small values degrade the statistics and large values miss too much of the fine-scale structure. We will use a value of ten with the smaller circle about 0.6% the size of the attractor and the larger circle about 6% the size of the attractor. Thus in Equation 3D we will use logarithms of base 10 (log10) instead of base 2 (log2).

Ideally the circles should be placed uniformly or randomly over the plane. However if we were to do this most of the circles would be empty and a very long calculation would be required to obtain an accurate estimate of the dimension. Instead we center the circles on the data points themselves. In this way the circles tend to enclose many points. However it represents a different type of average since it weighs more heavily the portions of the attractor where the points are most dense. Technically what we are calculating is called the correlation dimension since it involves the number of points that are correlated with each point in the data set. The correlation dimension is never greater than the fractal dimension but it tends not to be much smaller either.

The correlation dimension is only one of many ways to define the dimension of an attractor. The various methods differ in how the regions of the attractor are weighed in the average. It is probably the easiest method to implement and it gives more reliable results than the fractal dimension when the dimension of the attractor is greater than about two. The fractal dimension is also called the capacity dimension and it is closely related to the Hausdorff-Besicovitch dimension. Furthermore the correlation dimension is probably a more meaningful measure of the strangeness of the attractor since it includes information about its formation as well as its final appearance.

The number of data points required to provide an accurate estimate of the dimension is a question still being debated in the scientific literature. Therefore we will use an heuristic approach and continually update the dimension estimate with each iteration giving you an opportunity to decide when it seems to have converged to a unique value. To do this we must modify the procedure slightly. Rather than count the number of data points within a circle which would require that the calculation run to conclusion with the coordinates of all the points saved we will use the equivalent method of determining the probability that two randomly chosen points are within a certain distance of one another. To do this the distance of each new iterate from one of its randomly chosen predecessors is calculated. Now you see why we bothered to save the last 500 iterates! We will exclude the most recent 20 points since the iterates are likely to be abnormally highly correlated with their recent predecessors. Thus with each iteration we have only one additional calculation to do in which we compare the distance of the iterate to one of its randomly chosen predecessors and increment N1 and N2 as appropriate. The listing PROG07 shows the changes needed to calculate and display the fractal dimension.

PROG07 Changes Required in PROG06 to Calculate and Display the Fractal Dimension

1000 REM TWO-D MAP SEARCH (With Fractal Dimension)

1020 DIM XS(499)
YS(499)
A(504)
V(99)
1580 P% = 0: LSUM = 0: N = 0: NL = 0: N1 = 0: N2 = 0
1620 TWOD% = 2 ^ D%
2210 XS(P%) = X: YS(P%) = Y
2440 GOSUB 3900                 'Calculate fractal dimension
3030    LSUM = LSUM + LOG(DF): NL = NL + 1
3060    IF N > 1000 AND N MOD 10 = 0 THEN LOCATE 1
76: PRINT USING "##.##"; L;
3170 XL = XMIN - MX: XH = XMAX + MX: YL = YMIN - MY: YH = YMAX + 1.5 * MY
3190 YH = YH - .5 * MY
3400 LOCATE 1
1: PRINT CODE$
3420 LOCATE 1
63: PRINT "F =": LOCATE 1
73: PRINT "L ="
3900 REM Calculate fractal dimension
3910 IF N < 1000 THEN GOTO 4010         'Wait for transient to settle
3920 IF N = 1000 THEN D2MAX = (XMAX - XMIN) ^ 2 + (YMAX - YMIN) ^ 2
3930 J% = (P% + 1 + INT(480 * RND)) MOD 500
3940 DX = XNEW - XS(J%): DY = YNEW - YS(J%)
3950 D2 = DX * DX + DY * DY
3960 IF D2 < .001 * TWOD% * D2MAX THEN N2 = N2 + 1
3970 IF D2 > .00001 * TWOD% * D2MAX THEN GOTO 4010
3980    N1 = N1 + 1
3990    F = .434294 * LOG(N2 / (N1 - .5))
4000    LOCATE 1
66: PRINT USING "##.##"; F;
4010 RETURN

At this point you might wish to examine the fractal dimension of the various figures in this book as well as the dimension of those you create with PROG07. One thing you will notice is that the dimension of objects that resemble lines is often less than 1.0. One reason is that the points that make up the line are seldom uniformly distributed along its length. Remember that the correlation dimension is usually smaller than the fractal dimension. They are equal if the points are uniformly distributed over the attractor. The correlation dimension of a line consisting of a uniform distribution of points along its length would be exactly 1.0.

Also note that the dimension of most attractors varies considerably from one part of the attractor to another. Figure 3-11 is a good example of one in which parts of the attractor resemble thin lines and other parts resemble filled-in planes. It is obviously simplistic to characterize such an object by a single average dimension.

The dimension also depends on scale. It is properly defined in the limit where one zooms in very tightly on the attractor to observe its finest detail. However a calculation in this limit would take forever since an infinite number of iterations would be required to collect enough points to reveal the detail. Figure 3-13 is an example of an attractor that is nearly one-dimensional on a large scale but closer to two-dimensional on a fine scale. Our calculation provides what might be called a visual dimension since it is taken on a scale close to what the eye can visually resolve. In any case you should not ascribe undue significance to the calculated dimension.

Also note that we are using the word "dimension" to mean several different things. The maps that we are looking at are two-dimensional because they have two variables X and Y. However the attractor has a smaller dimension. We say the attractor is embedded in a two-dimensional space or that the embedding dimension is 2. A point or a line can be embedded in a plane but a ball cannot.

An attractor usually fills a negligible portion of the space in which it is embedded. That's why it's called an attractor! Points initially distributed throughout the embedding space are drawn to the attractor after a number of iterations and the remaining space is left empty. Thus the area of an attractor embedded in a two-dimensional space is zero and the volume of an attractor embedded in a three-dimensional space is zero and so forth.

It is also interesting that the fractal dimension and the Lyapunov exponents are not entirely independent. It has been conjectured but never completely proved that the fractal dimension is related to the two Lyapunov exponents by

F = 1 - L1 / L2 (Eq. 3E)

where L1 is the more positive of the two exponents and is the one we denote by L in the figures. If both Lyapunov exponents are known Equation 3E can be used to define a dimension of the attractor called the Lyapunov dimension. The Lyapunov dimension is also called the Kaplan-Yorke dimension after the scientists who made the conjecture.

This relation is reasonable since if the two exponents are equal but of opposite signs (L2 = - L1) the contraction in one direction is just offset by expansion in the other. A set of initial conditions spread out over a two-dimensional region thus maintains its area upon successive iteration. Such a mapping is said to be area-preserving or Hamiltonian after the nineteenth century Irish astronomer William Rowan Hamilton. On the other hand if the contraction is very rapid (L2 is large and negative) the initial conditions quickly collapse to a very elongated ellipse whose dimension is close to 1. Such a contraction is sometimes called filamentation.

Armed with information about the fractal dimension you can program the computer to be even more selective. For example the visually appealing attractors tend to have fractal dimensions slightly greater than 1 and thus you could reject those with smaller dimensions or those with dimensions close to 2. We will return to this intriguing possibility in Chapter 8.

3.5 Higher Order Disorder

With one-dimensional maps the attractors became more interesting when we considered terms higher than quadratic. It is straightforward to do the same with two-dimensional maps. For example the most general equations for two-dimensional cubic maps are:

Xn+1 = a1 + a2Xn + a3Xn2 + a4Xn3 + a5Xn2Yn

+ a6XnYn + a7XnYn2 + a8Yn + a9Yn2 + a10Yn3

(Eq. 3F)

Yn+1 = a11 + a12Xn + a13Xn2 + a14Xn3 + a15Xn2Yn

+ a16XnYn + a17XnYn2 + a18Yn + a19Yn2 + a20Yn3

Note that there are twenty coefficients which vastly increases the number of possible cases. The fourth-order case would have 30 coefficients and the fifth-order case would have 42 coefficients. If you prefer an equation a two-dimensional map of order O will have (O + 1)(O + 2) coefficients. We will code the cubic quartic and quintic cases with the letters F G and H respectively.

The changes that must be made to the program to generate attractors in two dimensions up to fifth order are given in the listing PROG08.

PROG08 Changes Required in PROG07 to Generate Attractors in Two Dimensions up to Fifth Order

1000 REM TWO-D MAP SEARCH (Polynomials up to 5th Order)

1020 DIM XS(499)
YS(499)
A(504)
V(99)
XY(4)
XN(4)
1060 OMAX% = 5                  'Maximum order of polynomial
1720 M% = 1: XY(1) = X: XY(2) = Y
1730 FOR I% = 1 TO D%
1740 XN(I%) = A(M%)
1750 M% = M% + 1
1760 FOR I1% = 1 TO D%
1770 XN(I%) = XN(I%) + A(M%) * XY(I1%)
1780 M% = M% + 1
1790 FOR I2% = I1% TO D%
1800 XN(I%) = XN(I%) + A(M%) * XY(I1%) * XY(I2%)
1810 M% = M% + 1
1820 IF O% = 2 THEN GOTO 1970
1830 FOR I3% = I2% TO D%
1840 XN(I%) = XN(I%) + A(M%) * XY(I1%) * XY(I2%) * XY(I3%)
1850 M% = M% + 1
1860 IF O% = 3 THEN GOTO 1960
1870 FOR I4% = I3% TO D%
1880 XN(I%) = XN(I%) + A(M%) * XY(I1%) * XY(I2%) * XY(I3%) * XY(I4%)
1890 M% = M% + 1
1900 IF O% = 4 THEN GOTO 1950
1910 FOR I5% = I4% TO D%
1920 XN(I%) = XN(I%) + A(M%) * XY(I1%) * XY(I2%) * XY(I3%) * XY(I4%) * XY(I5%)
1930 M% = M% + 1
1940 NEXT I5%
1950 NEXT I4%
1960 NEXT I3%
1970 NEXT I2%
1980 NEXT I1%
2000 NEXT I%
2010 M% = M% - 1: XNEW = XN(1): YNEW = XN(2)

The program PROG08 could have been written more compactly but it is done this way to simplify its extension to even higher dimensions. Examples of attractors produced by this program are shown in Figures 3-18 through 3-41.

FIRPGVTF

Figure 3-18. Two-Dimensional Cubic Map

FISHMQCH

Figure 3-19. Two-Dimensional Cubic Map

FJYCBMNF

Figure 3-20. Two-Dimensional Cubic Map

FLGROKJF

Figure 3-21. Two-Dimensional Cubic Map

FMGGNDPH

Figure 3-22. Two-Dimensional Cubic Map

FNHZBEET

Figure 3-23. Two-Dimensional Cubic Map

FNUYLCUR

Figure 3-24. Two-Dimensional Cubic Map

FOVFKWKE

Figure 3-25. Two-Dimensional Cubic Map

GFUXRRRU

Figure 3-26. Two-Dimensional Quartic Map

GGNXVYVA

Figure 3-27. Two-Dimensional Quartic Map

GLURFSRH

Figure 3-28. Two-Dimensional Quartic Map

GPFMQPPB

Figure 3-29. Two-Dimensional Quartic Map

GQDIDSBT

Figure 3-30. Two-Dimensional Quartic Map

GRMJQBCS

Figure 3-31. Two-Dimensional Quartic Map

GTPMJKFS

Figure 3-32. Two-Dimensional Quartic Map

GUETJGII

Figure 3-33. Two-Dimensional Quartic Map

HGEQGOYI

Figure 3-34. Two-Dimensional Quintic Map

HHVOIEGI

Figure 3-35. Two-Dimensional Quintic Map

HMSMTNCO

Figure 3-36. Two-Dimensional Quintic Map

HQBKSKIX

Figure 3-37. Two-Dimensional Quintic Map

HQDHFCND

Figure 3-38. Two-Dimensional Quintic Map

HSARYDPN

Figure 3-39. Two-Dimensional Quintic Map

HVHDXLMS

Figure 3-40. Two-Dimensional Quintic Map

HVNTBSGW

Figure 3-41. Two-Dimensional Quintic Map

Perhaps this is a good point to pause and reiterate in what sense these objects are attractors. If you choose initial values of X and Y somewhere near the attractor within its basin of attraction and substitute these values into the equations that describe the attractor the new values of X and Y will represent a point in the plane that is closer to the attractor. After a number of iterations the point will work its way to the attractor and thereafter it will move around on the attractor in some complicated manner eventually visiting every part of the attractor. The next position can always be simply and accurately predicted from the current position but the small inevitable uncertainty in position continually increases so that a long-term prediction is impossible except to say that the point will be somewhere on the attractor. You can think of the attractor as the set of all possible long-term solutions of the equations that produced it.

Besides the error in knowing perfectly the initial conditions there are also computer round-off errors at each iteration. Given the extreme sensitivity to small errors you may wonder whether any computer is capable of calculating correctly such a chaotic process. It is true that if the same chaotic equations are iterated on two computers using different precision or round-off methods the sequence of iterates will almost certainly be completely different after a few iterations. However the general appearance of the attractor will probably be the same. In such a case we say that the solution is structurally stable or robust. Furthermore according to the shadowing lemma an appropriate small change in initial conditions will produce a chaotic sequence that follows arbitrarily close to the computed one.

Since computers always round the results of calculations to a finite number of digits (or more precisely bits) there is a limited number of allowed values. Thus successive iteration of a map will always eventually repeat a previously obtained value whereupon the solution will reproduce exactly the same sequence of states as it did before. Strictly speaking every such solution is periodic and true chaos cannot be observed with a computer. However with double precision floating point variables which are normally 64 bits there are 264 or about 1019 possible values. Until the number of iterations approaches this value there is little cause to worry. For maps higher than one-dimensional this problem is even less serious since all the variables have to reach a previously existing state at the same time.

It is also interesting to realize that infinitely many periodic solutions are embedded in each attractor. These solutions are called periodic orbits. From wherever you start on the attractor you will eventually return to a point arbitrarily close to the starting point. This result is called the Poincaré recurrence theorem after Jules-Henri Poincaré a French mathematician who a hundred years ago portended the modern era of chaos. Thus by making only a small change in the starting point it is possible in principle to return exactly to the starting point which implies a periodic orbit with a period equal to the number of iterations required to return. Most of these orbits have very large periods however.

Every point on the attractor is arbitrarily close to such a periodic orbit but the chance that a randomly chosen point on the attractor lies on such an orbit is infinitesimal. We say that the periodic orbits are dense on the attractor. These orbits though infinite in number constitute a Cantor set of measure zero. The periodic orbits are unstable in the sense that if you get just slightly off the orbit you will continue to get farther away with each iteration.

The strange attractors exhibited in this book are examples of orbital fractals. They should be distinguished from escape-time fractals which show the basin of attraction and typically display with color the number of iterations required for points outside the basin to escape beyond some pre-defined region. The Mandelbrot and Julia sets are perhaps the best-known escape-time fractals. Escape-time fractals require much longer computing times to develop but provide dazzling displays with exotic fine-scale structure.

3.6 Strange Attractor Planets

The previous figures have obvious beauty but they generally lack symmetry. Nature mixes symmetry with disorder and our sense of beauty has developed accordingly. The Earth viewed from outer space is beautiful in part because the irregular features of the clouds and continents are superimposed on a nearly perfect sphere.

There are many ways to do the same with our attractors. Suppose for example X and Y are not the horizontal and vertical position in a plane but rather the longitude and latitude on the surface of the Earth. The result is an object that might resemble a strange planet with swirling clouds oceans canals craters and other features.

Note that mapping a plane onto a sphere is a nonlinear transformation. You can't wrap a piece of paper around a globe without a large non-uniform stretching. That's why Greenland looks larger than South America on most flat maps. When a sphere is projected onto a flat computer screen or onto the page of a book it is stretched so as to magnify the central portion of the attractor and compress the edges.

If q is the longitude (measured from zero at the right edge) and f is the latitude (measured from zero at the top) the X and Y coordinates of the projection of a sphere onto the screen are given by

Xp = cos q sin f

(Eq. 3G)

Yp = cos f

We get q from X by a scaling that keeps q in the range of 0 to p radians (180 degrees) since there is no need to plot points that lie on the back side of the planet. Similarly we get f from Y by a scaling that keeps f in the range of 0 (at the North Pole) to p radians (at the South Pole). The program modifications required to accomplish this transformation are given in the listing PROG09. The program allows you to toggle back and forth between the two types of projection by pressing the P key.

PROG09 Changes Required in PROG08 to Project Attractor onto a Sphere

1000 REM TWO-D MAP SEARCH (Projected onto a Sphere)

1110 PJT% = 1                   'Projection is spherical
2260    IF PJT% = 1 THEN GOSUB 4100     'Project onto a sphere
3200 XA = (XL + XH) / 2: YA = (YL + YH) / 2
3310 IF PJT% <> 1 THEN LINE (XL
YL)-(XH
YH)

B
3320 IF PJT% = 1 THEN CIRCLE (XA
YA)
.36 * (XH - XL)
3330 TT = 3.1416 / (XMAX - XMIN): PT = 3.1416 / (YMAX - YMIN)
3750 IF Q$ = "P" THEN PJT% = (PJT% + 1) MOD 2: T% = 3: IF N > 999 THEN N = 999
4100 REM Project onto a sphere
4110 TH = TT * (XMAX - XP)
4120 PH = PT * (YMAX - YP)
4130 XP = XA + .36 * (XH - XL) * COS(TH) * SIN(PH)
4140 YP = YA + .5 * (YH - YL) * COS(PH)
4150 RETURN

Figures 3-42 through 3-57 show some examples of two-dimensional attractors projected onto a sphere. Note that the features on the attractors tend to converge at the poles at the top and bottom of the figures. This convergence could be suppressed by using an area-preserving transformation that stretches the Y-values near the poles by the same factor that the X-values are compressed. This simplest way to produce this effect is just to delete line 4140.

ECSRKVVQ

Figure 3-42. Two-Dimensional Quadratic Map Projected onto a Sphere

ECVQKGHQ

Figure 3-43. Two-Dimensional Quadratic Map Projected onto a Sphere

EKPNERVO

Figure 3-44. Two-Dimensional Quadratic Map Projected onto a Sphere

EUWACXDQ

Figure 3-45. Two-Dimensional Quadratic Map Projected onto a Sphere

FKAWYMKA

Figure 3-46. Two-Dimensional Cubic Map Projected onto a Sphere

FLQOBBRS

Figure 3-47. Two-Dimensional Cubic Map Projected onto a Sphere

FLUCBPVB

Figure 3-48. Two-Dimensional Cubic Map Projected onto a Sphere

FMEGVTLM

Figure 3-49. Two-Dimensional Cubic Map Projected onto a Sphere

GJCQPYDV

Figure 3-50. Two-Dimensional Quartic Map Projected onto a Sphere

GLQGRLUF

Figure 3-51. Two-Dimensional Quartic Map Projected onto a Sphere

GNWCAVJO

Figure 3-52. Two-Dimensional Quartic Map Projected onto a Sphere

GTNSTDPQ

Figure 3-53. Two-Dimensional Quartic Map Projected onto a Sphere

HEAUYOII

Figure 3-54. Two-Dimensional Quintic Map Projected onto a Sphere

HFJFWKAS

Figure 3-55. Two-Dimensional Quintic Map Projected onto a Sphere

HLTQSSRK

Figure 3-56. Two-Dimensional Quintic Map Projected onto a Sphere

HOKEFWLH

Figure 3-57. Two-Dimensional Quintic Map Projected onto a Sphere

If you are using PowerBASIC or its predecessor TurboBASIC and VGA graphics you will note a slight incompatibility with the CIRCLE command that causes the size of the circle that surrounds the attractor to vary from case to case. In these dialects of BASIC the radius of the circle in SCREEN modes 11 and 12 is specified in units of the screen height rather than its width. If you encounter this problem try replacing the .36 * (XH - XL) in line 3320 with .5 * (YH - YL).

Planes and spheres are not the only two-dimensional surfaces onto which attractors can be projected. A cylinder is another possibility. The cylinder can be oriented with its axis either horizontal or vertical or tilted at some arbitrary angle. A torus is another possibility. You may be able to think of other more exotic surfaces onto which the attractors can be projected.

3.7 Designer Plaids

It is interesting that all the one-dimensional maps described in the previous chapter are included in the two-dimensional cases. One needs only to set the coefficients of the Y-equation to zero. For example a two-dimensional map equivalent to the logistic equation is given by the code EMu%M9. However since Y doesn't change with successive iterations a graph of Y versus X is simply a straight horizontal line.

To display the logistic parabola we need to replace X with the next iterate of X and Y with the second next iterate of X. Two successive iterations of a quadratic map requires a fourth-order equation. A code that will accomplish this is GMu%M13NHUIM10.

There are other examples of two-dimensional maps that are really one-dimensional maps in disguise. Suppose Xn+1 depends only on Yn and Yn+1 depends only on Xn. Then Xn+2 will depend only on Xn and we have a one-dimensional map for X in which Y is merely an intermediate value of X. The most general fifth-order polynomial example of such a case is

Xn+1 = a1 + a17Yn + a18Yn2 + a19Yn3 + a20Yn4 + a21Yn5

(Eq. 3H)

Yn+1 = a22 + a23Xn + a24Xn2 + a25Xn3 + a26Xn4 + a27Xn5

This case can be achieved by setting the remaining thirty coefficients to zero in PROG09 by adding the following line after line 2730:

2735 IF (I% > 1 AND I% < M% / 2 - O%) OR I% > M% / 2 + O% + 1 THEN MID$(CODE$ I% + 1 1) = "M"

The result is to produce a 25th-order one-dimensional polynomial map displayed in two dimensions.

Figures 3-58 through 3-62 show sample attractors obtained in this way. Notice that they fill in rectangular regions resembling a plaid tartan in sharp contrast to all the previous cases. These attractors are especially appropriate for projecting onto a sphere since the features line up east-west along parallels and north-south along meridians. Figures 3-63 and 3-64 show some examples of plaid planetary attractors.

ECMMMEWH

Figure 3-58. Two-Dimensional Quadratic Plaid Map

FNMMMMMM

Figure 3-59. Two-Dimensional Cubic Plaid Map

GEMMMMMM

Figure 3-60. Two-Dimensional Quartic Plaid Map

HKMMMMMM

Figure 3-61. Two-Dimensional Quintic Plaid Map

ERMMMQEA

Figure 3-62. Two-Dimensional Quadratic Plaid Map on a Sphere

HOMMMMMM

Figure 3-63. Two-Dimensional Quintic Plaid Map on a Sphere

You might wish to try adding colors to emulate a decorative cloth pattern. On way to do this is to color the pixels according to the number of times they are visited by the orbit. This is easily done by changing line 2300 in the program to

2300 PSET (XP YP) (POINT (XP YP) + 1) MOD 16

which causes the color of the existing point at (XP YP) to be tested and then plotted with the next color in the palette of 16 colors. In the next chapter we will discuss other ways to produce colorful attractors.

3.8 Strange Attractors that Don't

From the foregoing discussion you might conclude that all chaotic equations produce strange attractors. Such is not the case. Under certain conditions the successive iterates of an equation will wander chaotically throughout a region of the plane. There is no basin of attraction and initial conditions near but outside the chaotic region are not drawn to the region but rather lie on closed curves. Although the chaotic region is not a strange attractor it may have considerable beauty.

For a chaotic solution not to attract the area occupied by a cluster of nearby initial conditions must remain the same with successive iterations. The cluster will generally contract in one direction and expand in the other but the contraction and expansion just cancel producing a long thin filament of constant area. A characteristic of such a case is that the two Lyapunov exponents are equal in magnitude but of opposite signs. Such a system is said to be area-preserving or Hamiltonian.

You might think that Hamiltonian systems are relatively rare in nature since they require a special condition. However there are many important examples of Hamiltonian chaos. They arise because there are quantities in nature such as energy and angular momentum that in the absence of friction remain accurately constant no matter how complicated the behavior of the system. We say these quantities are conserved or that they are constants of the motion. The motion of a planet orbiting a binary-star system or the motion of an electron near the null in a magnetic field exhibit Hamiltonian chaos. A more familiar example is the filamentation of milk when it is stirred into coffee in which case the volume of the milk is conserved because liquids are nearly incompressible.

With equations such as those we have been using with randomly chosen coefficients the chance of inadvertently finding an area-preserving solution is essentially zero. However by placing appropriate conditions on the coefficients we can guarantee such solutions. An example of an area-preserving two-dimensional polynomial map is the following:

Xn+1 = a1 + a2Xn + a3Xn2 + a4Xn3 + a5Xn4 + a6Xn5 ± Yn

(Eq. 3I)

Yn+1 = a22 ± Xn

This map is fifth-order to provide seven arbitrary coefficients that ensure a large number of solutions. The coefficient labels are consistent with the general two-dimensional fifth-order map in which 33 of the coefficients have been set to zero. The two terms preceded with ± (a17 and a23) have coefficients of either +1 or -1 and this feature guarantees an area-preserving solution. If the signs are the same (both plus or both minus) chaotic solutions are not found. Hamiltonian chaos can occur when the signs are opposite. The negative product of these two coefficients is the Jacobian of the map (J = -a17a23). The Jacobian is a measure of the net contraction and it must equal 1.0 for a Hamiltonian system.

Hamiltonian cases can be produced by adding to PROG09 after line 2730 the following lines:

2735 IF I% > O% + 1 AND I% <> M% / 2 + 1 THEN MID$(CODE$ I% + 1 1) = "M"

2736 MID$(CODE$ M% / 2 - O% + 2 1) = "W": MID$(CODE$ M% / 2 + 3 1) = "C"

Sample Hamiltonian chaotic solutions are shown in Figures 3-64 through 3-71. Most of the cases resemble chains of islands in which each island contains a fixed point surrounded by closed contours that are not shown. These cases were produced using initial val