My Christmas Gift: Mandelbrot Set Computation In PythonMy mother likes fractals for their strange beauty. I decided to give her a simple way to generate beautiful fractal images of the Mandelbrot set. As you may guess, I'll be using Python for this. There are available Python code for this on the web, but they are either slow, or they don't produce nice images. Hence my own attempt at it. The code used here is available in a notebook on github or on nbviewer. I explore various ways to speed that code in How To Quickly Compute The Mandelbrot Set In Python . In order to tease you, here is one sample of what you'll be able to create very soon. It is an image of the Mandelbrot set. You can click on each image in this post to get a higher precision image. Mandelbrot SetYou can skip this explanation and directly go to the code if you wish, but I feel compelled to define precisely what we are computing here. Each point in the plane of coordinates (x,y) can be interpreted as a complex number with real part x and imaginary part y. Let c be such point. We define a series of complex numbers for it: z_{0} = 0 z_{1} = c z_{2} = z_{1}^{2} + c ... z_{n+1} = z_{n}^{2} + c
If the size z_{n} of these complex numbers stays bounded as n tends to infinity, then c belongs to the Mandelbrot set. Computing ImagesIn order to compute images, people usually compute the series until either z_{n} exceed a given horizon value, or n reaches a max number of iterations. In the later case, we assume that the point belongs to the Mandelbrot set. The following function computes that number with an horizon equal to 2. It can be shown that using any value greater or equal to 2 for the horizon yields the same set. import numpy as np from numba import jit The @jit says that I want the function to be compiled with Numba. See How The following function computes the number of iterations for all points in a given region of the plane. It gets the number of pixels to compute as arguments. @jit def mand It is rather fast. Computing a 1M pixels image takes about 1.5 seconds on my old laptop. DisplayWe will use matplotlib in a Jupyter notebook. We need to import it and we use the matplolib inline magic to display in the notebook. Most of the code is about decorating the image. The key piece is the use of the imshow function that displays the input array as an image. The function gets the coordinates of the region to compute plus the size of the image to compute, in inches. from matplotlib import pyplot as plt from matplotlib import colors %matplotlib inline
The last line is optional, only use it if you also want to save the image as a file. I used it to save the images shown in this entry with this function: image_counter = 30 def save_image(fig): global image_counter filename = "man It can be shown that the Mandelbrot set is entirely contained in the region where 2.5 <= x <= 0.5 and 1.25 <= y <= 1.25. Therefore, we call our function with these values:
It yields this image of the Mandelbrot set. Color MapsWhen I first shown this image to my wife, she asked if we could not use other colors. I must tell you that she is a painter, and that color choices is her bread and butter. Matplotlib lets us select the color palette when rendering an image. In order to produce the above image, matplotlib maps each value in the input array to a color, before rendering the result. The way values are mapped to colors is defined by color maps. A list of predefined color maps can be found here. We modify our function to accept a color map as input. Jet is the default color map. def mand My wife looked at all color maps, and she decided that the 'hot' color map was producing beautiful colors. mand yields The image is too dark and not all colors of the palette are used. Could it be that the pixel values aren't distributed correctly? There is only one way to answer that question: plot a histogram of the values taken by the pixels. The higher the bar in the histogram, the more pixels have this value.
yields We see that most pixels have a value between 0 and 10. We can spread the values using a low power normalization, for instance taking the 0.3 power of the values:
This looks better. We can achieve this normalization using a normalize object in our function. This object takes as argument the power coefficient. We pass it as the gamma argument to our function. def mand Executing mand
Much better, but the image isn't perfect. The region is segmented by color levels. Some may find it aesthetic, but I'd rather have smooth transitions between color levels. Anti AliasingWe can achieve smooth color transitions with the following function. It computes a floating point value for each pixel rather than an integer. The math behind it are given at the end of the post for those interested. The main difference with the previous function is that we use a very large horizon instead of 2, namely 2^{40} . @jit def mand Using it yields this image of the Mandelbrot set This is much better, we no longer see the sharp transitions we had between colors. Zooming InWe can now zoom in the image if we wish to, by modifying the coordinates of the region we want to display. For instance, this zooms into the valley between the two largest black areas: mand Oops, the result does not look right. We need to use a larger number of iterations for a more accurate result, say 1024 iterations. Win. For the images below we use 2048 iterations. Computation time stays reasonable, a couple of seconds for a 1M pixels image, to my surprise. In fact, the less black in the picture, the faster the computation. I'll let readers figure out why. Zooming on the right side of the valley yields a nice set of inverted sea horses. Zooming on one of these sea horses Zooming on its tail Zooming on the tiny black dot near the lower left corner Surprise! The black dot appears to be a copy of the whole Mandelbrot set itself. This is not by chance. It is precisely because small parts are similar to the whole that the Mandelbrot set is a fractal. Your TurnYou are all set. You can use the code to zoom in any region you'd like. Beware though that there will be a limit on how far you can go before hitting rounding error problems. You can also explore using other color maps. For instance, here is the Mandelbrot set with the magma color map. And here is the latest detail with the gnuplot2 color map. You may need to change the gamma parameter to get better results, depending on the color map and on the number of iterations. You can also create larger images, by providing new values for the width and height. These are expressed in inches. For instance, to create a 20x20 inches image you can call: mand If you want to save images for a printer, then you should change the dpi value in the function from 72 to 300. Beware, this will slow computation. As promised I'll conclude with an explanation of why our second mandelbrot_set function provides smooth colors. Non mathematically oriented readers can skip it safely. Some MathLet's see how we moved to an integer number of iterations to a floating point value. The first step is to chose a very large horizon h. We use h = 2^{40} which is about one trillion. For a given point c, we have z_{n+1} = z_{n}^{2} + c When the size of z_{n} is approaching h, then we can neglect c, and we have that z_{n+1} = z_{n}^{2} Taking the log of the size of on both side we get log( z_{n+1} ) = 2 log( z_{n} ) Therefore, log( z_{n} ) grows like the power of 2: log( z_{n}  ) = 2^{n+a} for some constant a. Taking the log base 2 we get f(n) = n+a where f(n) = log( log( z_{n} ) ) / log(2) Let's define g(x) = x + a for any real number x. We have that g(x) = f(x) when x is integer. Let N be the first iteration where z_{N} exceeds the horizon h: z_{N1} < h <=  z_{N}  Taking the log twice g(N1) < log( log(h) ) / log(2) <= g(N) Instead of returning N, we want to return A such that g(A) = log( log(h) ) / log(2) A is the value where g crosses the horizon. We would then have g(N1) < g(A)<= g(N) which is what we are looking for. We have A = g(A)  a = log( log(h) ) / log(2)  a = log( log(h) ) / log(2)  f(N) + N This is what our function computes. Q.E.D. I didn't invent that approach, but my exposure may be different from others. Wikipedia provides these references:
