Plotting in Python

Problems taken from Computational Physics by Newman Chapter 3

2. Curve plotting

Although the \verb|plot| function is designed primarily for plotting standard xy graphs, it can be adapted for other kinds of plotting as well.
Make a plot of the so-called deltoid curve, which is defined parametrically by the equations

x = 2\cos \theta + \cos 2 \theta
y = 2 \sin \theta - \sin 2 \theta

where 0\le\theta<2\pi. Take a set of values of \theta between zero and 2\pi and calculate x and y for each from the equations above, then plot y as a function of x. Taking this approach a step further, one can make a polar plot r=f(\theta) for some function f by calculating r for a range of values of \theta and then converting r and \theta to Cartesian coordinates using the standard equations x = r\cos\theta, y = r\sin\theta. Use this method to make a plot of the Galilean spiral r = \theta^2 for 0\le\theta\le10\pi.
Using the same method, make a polar plot of “Fey’s function”

r = e^{\cos\theta} - 2 \cos 4\theta + \sin^5 \frac{\theta}{12}

in the range 0 \le \theta \le 24\pi.

3.6 Deterministic chaos and the Feigenbaum plot

One of the most famous examples of the phenomenon of chaos is the logistic map, defined by the equation

x' = rx(1-x).

For a given value of the constant r you take a value of x – say x=\frac{1}{2} – and you feed it into the right-hand side of this equation, which gives you a value of x’. Then you take that value and feed it back in on the right-hand side again, which gives you another value, and so forth. This is a iterative map. You keep doing the same operation over and over on your value of x, and one of three things happens:

  • The value settles down to a fixed number and stays there. This is called a fixed point. For instance, x=0 is always a fixed point of the logistic map. (You put x=0 on the right-hand side and you get x’=0 on the left.)
  • It doesn’t settle down to a single value, but it settles down into a periodic pattern, rotating around a set of values, such as say four values, repeating them in sequence over and over. This is called a limit cycle.
  • It goes crazy. It generates a seemingly random sequence of numbers that appear to have no rhyme or reason to them at all. This is deterministic chaos. “Chaos” because it really does look chaotic, and “deterministic” because even though the values look random, they’re not. They’re clearly entirely predictable, because they are given to you by one simple equation. The behavior is determined, although it may not look like it.

Write a program that calculates and displays the behavior of the logistic map. Here’s what you need to do. For a given value of r, start with x=\frac{1}{2}, and iterate the logistic map equation a thousand times. That will give it a chance to settle down to a fixed point or limit cycle if it’s going to. Then run for another thousand iterations and plot the points (r,x) on a graph where the horizontal axis is r and the vertical axis is x. You can either use the plot function with the options “ko” or “k.” to draw a graph with dots, one for each point, of you can use the scatter function to draw a scatter plot (which always uses dots). Repeat the whole calculation for values of r from 1 to 4 in steps of 0.01, plotting the dots for all values of r on the same figure and then finally using the function show once to display the complete figure.

Your program should generate a distinctive plot that looks like a tree bent over onto its side. This famous picture is called the Feigenbaum plot, after its discoverer Mitchell Feigenbaum, or sometimes the figtree plot, a play on the fact that it looks like a tree and Feigenbaum means “figtree” in German.

Give answers to the following questions:
For a given value of r what would a fixed point look like on the Feigenbaum plot? How about a limit cycle? And what would chaos look like? Based on your plot, at what value of r does the system move from orderly behavior (fixed points or limit cycles) to chaotic behavior? This point is sometimes called the “edge of chaos.”

The logistic map is a very simple mathematical system, but deterministic chaos is seen in many more complex physical systems also, including especially fluid dynamics and the weather. Because of its apparently random nature, the behavior of chaotic systems is difficult to predict and strongly affected by small perturbations in outside conditions. You’ve probably heard of the classic exemplar of chaos in weather systems, the butterfly effect, which was popularized by physicist Edward Lorenz in 1972 when he gave a lecture to the American Association for the Advancement of Science entitled, “Does the flap of a butterfly’s wings in Brazil set off a tornado in Texas?” (Although arguably the first person to suggest the butterfly effect was not a physicist at all, but the science fiction writer Ray Bradbury in his famous 1952 short story A Sound of Thunder, in which a time traveler’s careless destruction of a butterfly during a tourist trip to the Jurassic era changes the course of history.

Comment: There is another approach for computing the Feigenbaum plot, which is neater and faster, making use of Python’s ability to perform arithmetic with entire arrays. You could create an array r with one element containing each distinct value of r you want to investigate: [1.0, 1.01, 1.02, … ]. Then create another array x of the same size to hold the corresponding values of x, which should all be initially set to 0.5. Then an iteration of the logistic map can be performed for all values of r at once with a statement of the form x = r*x*(1-x). Because of the speed with which Python can perform calculations on arrays, this method should be significantly faster than the more basic method above.

3.7 The Mandelbrot set

The Mandelbrot set, named after its discoverer, the French mathematician Benoit Mandelbrot, is a fractal, an infinitely ramified mathematical object that contains structure within structure within structure, as deep as we care to look. The definition of the Mandelbrot set is in terms of complex numbers as follows. Consider the equation

z' = z^2 + c,

where z is a complex number and c is a complex constant. For any given value of c this equation turns an input number z into an output number z’. The definition of the Mandelbrot set involves the repeated iteration of this equation: we take an initial starting value of z and feed it into the equation to get a new value z’. Then we take that value and feed it in again to get another value, and so forth. The Mandelbrot set is the set of points in the complex plane that satisfies the following definition:

For a given complex value of c, start with z=0 and iterate repeatedly. If the magnitude |z| of the resulting value is ever greater than 2, then the point in the complex plane at position c is not in the Mandelbrot set, otherwise it is in the set.

In order to use this definition one would, in principle, have to iterate infinitely many times to prove that a point is in the Mandelbrot set, since a point is in the set only if the iteration never passes |z|=2 ever. In practice, however, one usually just performs some large number of iterations, say 100, and if |z| hasn’t exceeded 2 by that point then we call that good enough.

Write a program to make an image of the Mandelbrot set by performing the iteration for all values of c=x+ iy on an N\times N grid spanning the region where -2\le x\le 2 and -2\le y\le 2. Make a density plot in which grid points inside the Mandelbrot set are colored black and those outside are colored white. The Mandelbrot set has a very distinctive shape that looks something like a beetle with a long snout—you’ll know it when you see it.

Hint: You will probably find it useful to start off with quite a coarse grid, i.e., with a small value of N — perhaps N=100 — so that your program runs quickly while you are testing it. Once you are sure it is working correctly, increase the value of N to produce a final high-quality image of the shape of the set.

If you are feeling enthusiastic, here is another variant of the same exercise that can produce amazing looking pictures. Instead of coloring points just black or white, color points according to the number of iterations of the equation before |z| becomes greater than 2 (or the maximum number of iterations if |z| never becomes greater than 2). If you use one of the more colorful color schemes Python provides for density plots, such as the “hot” or “jet” schemes, you can make some spectacular images this way. Another interesting variant is to color according to the logarithm of the number of iterations, which helps reveal some of the finer structure outside the set.