More with the Mandelbrot Set

The Mandelbrot Set is the set of points for which iterating the equation

  zn+1=zn2+z0

does not diverge to infinity, but rather stays finite.

An obvious question arises: how does it behave as it stays finite?

A fixed point

It could simply converge to a finite value. If so, that value should be easy to find:

z=z2+z0
z=0.5±√(0.25-z0)

Such a value is called a fixed point of the equation, for, once one reaches the point, one remains fixed there. But will this fixed point attract towards itself, or repel? In other words, if the fixed point is called zf, does a close-by point, zf+ε move closer, or further away, as iterations progress?

If M(z)=z2+z0, the question can be re-expressed as whether

|M(zf+ε)-M(zf)|<ε

But |M(z+ε)-M(z)| is simply εdM/dz, and dM/dz is simply 2z.

So the condition for an attractive fixed point to exist is

|1±√(1-4z0)|<1

This is most easily solved by considering the boundary where the modulus equals one. If a (complex) number's modulus is one, then it can be written as exp(it), so we can solve

1±√(1-4z0)=exp(it)
z=0.5exp(it)(1-0.5exp(it))

Of course one might not recognise that parametric equation. So python to the rescue!

#!/usr/bin/python3
import numpy as np
import matplotlib.pyplot as plt

t=np.linspace(0,2*np.pi,100)

z=0.5*np.exp(1j*t)*(1-0.5*np.exp(1j*t))
x=z.real
y=z.imag

plt.plot(x,y)
plt.axis('equal')
plt.show()

cardioid

It is a cardioid. A more usual form, (try proving that these are equivalent, with z=x+iy) might be:

  x=0.25+0.5cos(φ))[1-cos(φ)]
  y=0.5sin(φ)[1-cos(φ)]

And the area of this part is 3π/8 (c. 1.178). How does this compare with any area calculated following the suggestion on the Mandelbrot page? The cusp of a cardioid is usually located at the origin, but the extra constant term in the expression for x moves it to x=0.25.

Note that this is the condition that an attractive fixed point exists. It has not been shown that the sequence starting at zero will converge to it, rather than the attractive fixed point at infinity. In practice it does.

Cycles

Another possibility is that the sequence converges towards some form of cycle. Consider z0=-1. There is now a fixed cycle of -1, 0, -1, 0, ..., a cycle of period two. Other values of z0 close to -1 also behave in this way. Indeed, all points within a circle of radius 0.25 centred on -1 end up in an attractive cycle of period two. The area of this part is, of course, π/16.

One could write some python to iterate z2+z0 a few hundred times, and then print out the next twenty values to see if one can spot other cycle lengths. What happens around the following points?

#!/usr/bin/python3

z0=eval(input("Input z0 in the form a+bj: "))

z=0
for i in range(0,200):
    z=z*z+z0

for i in range(0,20):
    z=z*z+z0
    print(z)

The Extent of the Set

The Mandelbrot Set is often drawn with the real axis extending from -2 to +1, and the imaginary axis from -1.5i to 1.5i, resulting in a square window. One can fairly quickly demonstrate that, on the real line, the set extends from -2 to 0.25. But can you find numerically the how large a rectangle is required to contain the whole set?

(The greatest extent in the imaginary direction is known to be given by two of the roots of

x3+2x2+2x+2=0

but I am unaware of an analytic result for the greatest extent in the positive real direction, save that the cardioid above includes the points 0.375±√3i/8 and thus gives a lower bound.)

The Point at -0.75

The point at -0.75 lies on the boundary between the main cardioid in which iterations are attracted to a single fixed value, and a circle centred on -1 in which iterations are attracted to a cycle of period two. But how thick is this "neck" which joins the "head" on the Mandelbrot "beetle" to its "body"? Consider the following code for investigating the number of iterations for a starting value of -0.75+εi to diverge.

#!/usr/bin/python3

epsilon=0.1
print("epsilon   iters  epsilon*iters")

while (epsilon>=1e-6):
    z0=complex(-0.75,epsilon)
    i=1
    z=z0
    while (z.real**2+z.imag**2<4):
        z=z*z+z0
        i=i+1
    print("%7g  %7d  %7f"%(epsilon,i,i*epsilon))
    epsilon=0.1*epsilon

So the neck is infinitely thin at this point. This is quite a good way of checking how accurate a depiction of the Mandelbrot Set is. To improve the quality, one must increase the maximum number of iterations before one assumes that a sequence does not diverge to infinity.

The surprising constant product of ε and the number of iterations before divergence is clear was discovered computationally by Dave Boll in 1991, then a graduate student at Colorado State University studying Computer Science. The result was later proved analytically by Aaron Klebanoff (Fractals 9 393 (2001), also found here).

Why is >2 equivalent to infinity?

The test for a point diverging to infinity is that its modulus exceeds two. At this point the modulus of z2 will be more than two greater than the modulus of z. If one assumes that the modulus of z0 does not exceed two, then, even if z2 and z0 lie in opposite directions in the complex plane, it must be the case that |z2+z0|>|z|. So the modulus will increase on every iteration.

Ad infinitum?

Mathematics tells us that the boundary of the Mandelbrot Set never becomes smooth no matter how deeply one zooms in. Computers using standard double-precision arithmetic have a finite precision, and eventually one runs out of zoom. Consider the gallery image entitled "Surrounded by Spirals, off-centre detail zoomed." That is starting to reach the limits of double-precision arthmetic. If one zooms in on the centre one sees

full zoom

This is starting to look a little rough at the edges, and to show that we have run out of precision, consider a further zoom into the large "bud" at the top of the cardioid

over zoom

But why are the blocks which are appearing rectangular? Each block will represent those pixels which get rounded to the same double precision number. Computers store numbers in binary, with 53 digits. The x co-ordinates are around 0.8, and 0.8 in binary starts 0.11 (i.e. 1/2 + 1/4). The y cordinates are around 0.2, and 0.2 in binary starts 0.0011 (i.e. 1/8 + 1/16). With 53 significant digits, the x coordinates are of the form:

0.11xx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx x

but the y coordinates are of the form 

0.0011 xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxx

So for the coordinates chosen, the minimum value of ε for which y+ε can be distinguished from y is four times smaller than the minimum value for which x+ε can be distinguished from x, and so the blocks are rectangles which are four times as wide as they are tall.

If the same area is calculated using a code which goes beyond double precision, one can obtain

over zoom recalculated at high precision

This was produced using a fixed point code. The parameters were

(-0.814276764228418107+0.1901980148050127445i) to (-0.814276764228412794+0.1901980148050180577i)
maxiter=90000
bshift=80

and the code used was the last from the page of C-enhanced fixed-point Mandelbrot codes. Emulating configurable-precision fixed-point arithmetic has a huge performance cost over using the fixed-precision hardware floating point functionality of modern CPUs. These calculations were done on an Intel CPU of the Haswell generation, and the performance penalty for the fixed point version was a factor of fourteen.

A further example can be found on the 128 bit fixed point page.

sci-pi.org.uk Legal / Contact