Nostalgi-O-Vision, activate!

A month or so after I was born, my parents bought an Atari 400 game console. It plugged into the television set, and it had a keyboard with no moving keys, intended to be child- and spill-proof. Thanks to the box of cartridges we had beside it, Asteroids and Centipede were burnt into my brain at a fundamental level. The hours I lost blowing up all my own bases in Star Raiders — for which accomplishment the game awarded you the new rank of “garbage scow captain” — I hesitate to reckon. We also had a Basic XL cartridge and an SIO cassette deck, so you could punch in a few TV screens’ worth of code to make, say, the light-cycle game from TRON, and then save your work to an audio cassette tape.

From my vantage point in the twenty-first century, it seems so strange: you could push in a cartridge, close the little door, turn on your TV set and be able to program.

Granted, you had to work with a cumbersome programming language whose use might leave you with bad programming habits, but it was so immediate. Nowadays, we have more capable languages, ones which you can begin with as a novice and do serious work in later. (Ones in which you don’t have to put a number at the front of each line.) But I can’t just plug in a Python cartridge and get it on my TV; on many computers, I have to go download it and then fuss through the installation process, not only for the language itself but also for additional packages needed for particular applications. Said process is, of course, different, not just for each operating system, but for each incremental version thereof. I had to go through this for a roomful of people and their laptops several weekends ago, getting SciPy and matplotlib and NetworkX running. It wasn’t pretty, although the end result is quite useful. (And then we have to keep it going. Here in the office, the very next week, we saw all our SciPy/matplotlib/Sage installations break because, as far as we can tell, Mac OS X decided to upgrade itself and thereby clear-cut the whole dependency tree. Software only ever “just works” if you never do anything with it.)

But anyway! In the spirit of those long-forgotten BASIC code listings in algebra books, let’s have a bit of code with pædagogical utility. We’ve already seen how to implement Conway’s Game of Life in a few minutes; how about we try for another classic, the logistic map?

This is a mathematical idealization of a population of organisms, living in an environment which has a carrying capacity: the environment can’t sustain more than some fixed number of individuals. Here, we work in discrete time. That is to say, the population jumps from one size to another according to the rule we define, rather than changing smoothly; no time exists in between the ticks of the clock. We’ll let $x_t$ denote the fraction of the available space which is filled at time $t$. We specify the behaviour of our system by writing the update rule which tells us what $x_{t+1}$ will be, given $x_t$. In this case, we use

$$ x_{t+1} = rx_t(1 – x_t). $$

When $x_t$ is small, then this roughly approximated by ordinary exponential growth, $x_{t+1} = rx_t$, but if $x_t$ is closer to 1, then the second term within the parentheses kicks in, and the growth rate is diminished.

An aside: why is this called “logistic” growth? Nobody knows for sure why Verhulst first decided to call his equation of this sort “la courbe logistique“. He taught at Belgium’s military academy, so he may have picked up the term logistics, which military folks had just recently started to use in the modern sense, referring to how one goes about provisioning armed forces with food and other necessary materiel. This may have suggested logistique as a description for an equation which treats population growth in the case of limited resources. Also, in Verhulst’s day (1840s), logistique referred to mathematical work done through explicit arithmetical computations — what we might call “number-crunching” — so he may have chosen that term to stress the numerical emphasis of his approach to demography.

Anyway! The other day, I hacked together logistic.py, which implements the logistic map and plots the results. I wanted something which demonstrated (for the weekend class I mentioned earlier) a few features of the Python language which I’ve found useful: reading arguments in from the command line, using try/except blocks to handle unexpected circumstances and so forth. The whole thing took less time than I’ve spent figuring out what to say in prefacing it. And that includes making the pretty pictures.

First, let’s look at what the program makes by default. Our Python code follows a standard scheme for exploring the logistic map: pick a value of the growth rate $r$, pick some starting value of the population size $x_0$ and run through the update rule a bunch of times. This is to give the system a chance to settle down. Then, iterate again several times, plotting the population size. If the population size stops bouncing around and settles to a steady, fixed value, then all the dots we plot will lie on top of each other. If the population size bounces between two values, back and forth, cyclically, then we’ll see two dots.

The big picture:

When the growth-rate parameter $r$ is less than 1, the population size drops to zero and stays there. For larger values of $r$, we still have a fixed point, but its value increases with $r$. Then, weirdly, we get a cyclic behaviour pattern: $x$ jumps back and forth between two numbers. The separation between these two numbers grows as we crank up $r$. From two, the cycle then splits into four.

Soon enough, the behaviour of $x$ becomes a chaotic mess!

(The SciPy software provides different options for plotting points, so you can tweak the code to make thinner lines, different colours, etc. My choices in this demo are biased by what I think looks cool.)

More weirdly still, if we zoom into one of those branching points, we get the entire original picture, reproduced in miniature. SciPy offers a zoom feature in its plotting window; by passing different command-line arguments to logistic.py, we can tell it to draw a sub-section of the bifurcation diagram with the same amount of detail we used in the original.

Notice the axes on this plot:

We can continue the process. Zoom again!

The tiny piece of the diagram in between $r = 3.560$ and $r = 3.575$ reproduces the original, in shrunken form. This self-similarity is a defining characteristic of fractal shapes.

Homework: find manually, and/or code a program to determine, the ratios between the values of $r$ at the successive branch points.

(Footnote: yes: somebody did implement the logistic map in BASIC. Some people juggle geese . . .)