The two latest posts have been about the Mandelbrot set, the set of complex numbers c such that iterations of

f(z) = z² + c

remain bounded. It’s easy to see that the sequence of iterates will go off to infinity if at any step |z| > 2.

For each c, we can look at the escape time, the number of iterations it takes before an value has magnitude greater than 2. The Mandelbrot set is the set of points with infinite escape time.

In the plot below, I sampled 1,000,000 points at random and colored each according to the corresponding escape time. The lighter blue color corresponds to earlier escape.

The white points in the interior have an escape time of more than 10 (and possibly infinity).

When plotting the Mandelbrot set, how do you know whether a point will never escape? You don’t, not with certainty. But by looking at the distribution of escape times you can see that most points that are going to escape do so quickly. For example, after sampling 1,000,000 points, doing 2000 iterations for each point, 99.9% of points that escape do so within 248 iterations.

If we were to plot a histogram of the escape times it would simply look like a spike on the left end. Using a logarithmic scale for the vertical axis lets us see more of what’s going on. Note that the escape times decrease rapidly, even on a log scale.

By counting what portion of points do not escape, we can calculate a Monte Carlo estimate of the area of the Mandelbrot set, and in fact this is the best way to calculate the area. There is an exact expression for the area in terms of an infinite sum, but the sum converges so slowly that it is unusable; Monte Carlo sampling is faster.

Using 1,000,000 random values of c, a total of 905,444 points from the square [−2, 2] × [−2, 2] escape within 2,000 iterations. Since the square has area 16, this gives us an estimate

A = 16 × (1000000 − 905,444)/1000000 = 1.5129.

A 95% confidence interval for the area is

A ± 16 × 1.96 × (p (1 − p)/1000000) = (1.5082, 1.5176)

where p = 905444/1000000.

A study [1] found the area to be 1.506484 ± 0.000004, which is just outside of our confidence interval. Some of the difference is due to chance, and some of it is due to our stopping at 2000 iterations. If we had gone done more iterations, more points would escape (though not very many, for reasons explained above).

Update: I redid my experiment with 10,000,000 samples and 32,000 iterations and got an area estimate of 1.507954 with a 95% confidence interval (1.505056, 1.510851), which would include the estimate cited above. However, for reasons I give here I suspect the confidence interval in [1] is not as small as reported.

[1] See OEIS A098403

As mentioned in the previous post, most of the area in the Mandelbrot set comes from two regions. The largest is the blue cardioid region below and the next largest is the orange disk.

The blue cardioid is the set of points c such that iterations of z² + c converge to a fixed point. The orange disk is the set of points that converge to a cycle of period 2.

We can watch the convergence by picking a point and drawing a line between the iterates with the following Python function.

import matplotlib.pyplot as plt

def plot_iterates(c):
    f = lambda z: z**2 + c
    z = c
    s = [z]
    for i in range(1, 20):
        z = f(z)
        s.append(z)
        plt.plot([s[i-1].real, s[i].real], [s[i-1].imag, s[i].imag], color="C0")
    plt.gca().set_aspect("equal")
    plt.show()

First, we start with a point in the cardioid that converges to a fixed point. We could say it converges to a cycle of period 1.

plot_iterates(0.05 + 0.3j)

Next, lets start at a point in the orange disk with period 2.

plot_iterates(-1 + 0.2j)

We can find points with any period n. Here’s one of period 3.

plot_iterates(-0.13 + 0.66j)

And here’s one of period 4.

plot_iterates(0.256 + 0.55j)

The values of c that converge to a cycle with period q are “bulbs” tangent to the cardioid. There are φ(q) bulbs where φ is Euler’s totient function, i.e. φ(q) is the number of positive integers less than and relatively prime to q. The orange disk is the first bulb, the only one of its kind because φ(2) = 1. There are two bulbs with period 3 because φ(3) = 2.

For every p relatively prime to q there is a bulb tangent to the cardioid at w(1 − w) where

w = ½ exp(2πip/q)

and the diameter of the bulb is approximately 1/q.

These bulbs don’t account for everything in the Mandelbrot set. There’s a lot more to explore.

By the way, here is a plot of the Mandelbrot set created by selecting values of c at random and testing whether iterations of z² + c remain bounded.