by Anand | tags: [ fractal mandelbrot julia python tutorial mathematics visualization matplotlib numpy ]
Fractals with Python - The Mandelbrot Set
In parts I and II of the fractal posts, I wrote about the Julia fractal set and its concurrent implementation. In the third post of this series, I am going to talk about its close cousin, the Mandelbrot Set.
Recap - Julia Set
Let us recall the definition of the Julia set. It is a mathematical set of points $z_i$ in the complex plane which satisfies the condition,
$$z_{i+1} = { z_i^2 + c }$$
where $$z_i ≠ ∞$$
We implemented this using the following Python code.
import numpy as np
def julia_set(width, height, zoom=1, niter=256):
""" A julia set of geometry (width x height) and iterations 'niter' """
w,h = width, height
pixels = np.arange(w*h, dtype=np.uint16).reshape(h, w)
# Pick some defaults for the real and imaginary constants
# This determines the shape of the Julia set.
c_real, c_imag = -0.7, 0.27
for x in range(w):
for y in range(h):
# calculate the initial real and imaginary part of z,
# based on the pixel location and zoom and position values
zx = 1.5*(x - w/2)/(0.5*zoom*w)
zy = 1.0*(y - h/2)/(0.5*zoom*h)
for i in range(niter):
radius_sqr = zx*zx + zy*zy
# Iterate till the point is outside
# the circle with radius 2.
if radius_sqr > 4: break
# Calculate new positions
zy,zx = 2.0*zx*zy + c_imag, zx*zx - zy*zy + c_real
color = (i >> 21) + (i >> 10) + i * 8
pixels[y,x] = color
return pixels
The Mandelbrot Set
The Mandelbrot Set - named after the French mathematician Benoit Mandelbrot who studied its properties extensively, is intimately connected to the Julia set.
It is defined as the set of complex points c, for which the Julia set is connected. In simpler words, this means the set of points where the Julia set contains the origin (0, 0)
.
The Mandelbrot set then forms the set of co-efficients c
of the Julia set, which satisfies this criterion.
Before we implement the Mandelbrot set in Python, let us reimplement the Julia set using complex numbers.
Julia Set - Using complex numbers
Since Python has built-in support for complex numbers, the original Julia Set code can be reimplemented to take advantage of this to be simpler and visually closer to the original mathematical definition of the set. This is the code for the same.
import numpy as np
# Julia Set reimplemented using complex numbers
def julia_set(width, height, cx=-0.7, cy=0.27, zoom=1, niter=256):
""" A julia set of geometry (width x height) and iterations 'niter' """
w,h = width, height
pixels = np.arange(w*h, dtype=np.uint16).reshape(h, w)
# Pick some defaults for the real and imaginary constants
# This determines the shape of the Julia set.
c = complex(cx, cy)
for x in range(w):
for y in range(h):
# calculate the initial real and imaginary part of z,
# based on the pixel location and zoom and position values
zx = 1.5*(x - w/2)/(0.5*zoom*w)
zy = 1.0*(y - h/2)/(0.5*zoom*h)
z = complex(zx, zy)
for i in range(niter):
radius_sqr = abs(z)
# Iterate till the point is outside
# the circle with radius 2.
if radius_sqr > 4: break
# Calculate new positions
z = z**2 + c
color = (i >> 21) + (i >> 10) + i * 8
pixels[y,x] = color
return pixels
The changes should be fairly self explanatory to anyone who remembers their complex number algebra and understands the original mathematical definition of the Julia Set.
With this re-implementation, we can now derive the code for the Mandelbrot Set from it.
Mandelbrot Set Generator
The Mandelbrot Set as defined above, is the set of co-efficients c
for which the corresponding Julia sets are centered on the origin. The logic of calculating the next point remains the same.
Hence by definition, for every point (x, y)
in the real plane, we need to find the set of co-efficients (cx, cy)
for the complex co-efficient c
in the complex plane which satisfies this criteria.
Let us walk through the Julia Set code step by step and build the corresponding Mandelbrot Set code.
We need to do the imports just as before.
"""
mandelbrot.py - Implement the Mandelbrot set in Python using matplotlib & numpy,
by reusing code from the Julia set.
"""
import matplotlib.pyplot as plt
import numpy as np
The function definition has a slight change as it no longer needs to accept any parameters for the complex constant c
, as in this case c
always defaults to the origin. However we will use a similar numpy array to hold the pixel data so that part of the code is the same as before.
Note that we also pass two additional arguments namely x_off
and y_off
for the x and y offsets for the viewport. This is useful when you are zooming the fractal and want to pan to the areas you want to zoom.
def mandelbrot_set(width, height, zoom=1, x_off=0, y_off=0, niter=256):
""" A mandelbrot set of geometry (width x height) and iterations 'niter' """
w,h = width, height
pixels = np.arange(w*h, dtype=np.uint16).reshape(h, w)
Now on to the main logic, the for loop that iterates over the dimensions (w x h)
.
# The mandelbrot set represents every complex point "c" for which
# the Julia set is connected or every julia set that contains
# the origin (0, 0). Hence we always start with c at the origin.
for x in range(w):
for y in range(h):
# calculate the initial real and imaginary part of z,
# based on the pixel location and zoom and position values
# We use (x-3*w/4) instead of (x-w/2) to fully visualize
# the fractal along the x-axis
zx = 1.5*(x + x_off - 3*w/4)/(0.5*zoom*w)
zy = 1.0*(y + y_off - h/2)/(0.5*zoom*h)
z = complex(zx, zy)
c = complex(0, 0)
The main difference from the Julia Set in this section is that we initialize the co-efficient c
to (0, 0)
(the origin) every time inside the loop, as per the definition.
Now onto the inner loop that performs the iterations.
for i in range(niter):
if abs(c) > 4: break
# Iterate till the point is outside
# the circle with radius 2.
# Calculate new positions
c = c**2 + z
In the code for Julia set we iterated over the point z
as $ z = z^2 +c$ and stopped when it fell out of the circle with radius 2.
Here since we are interested in the co-efficient space, we need to iterate over that. In other words, the roles of z
and c
are switched for the Mandelbrot Set, when compared to the Julia Set.
The rest of the code is the same as before, with the difference of the color palette. We want to enhance the blues in this color palette so the code is a bit different.
color = (i << 21) + (i << 10) + i * 8
pixels[y,x] = color
return pixels
Here is the full code, including the function which displays the generated pixels using matplotlib
. (Note that we have turned off the display of the co-ordinate axes and also accept a cmap parameter - for the colormap of the plot - in this version of the code.)
import matplotlib.pyplot as plt
import numpy as np
def mandelbrot_set(width, height, zoom=1, x_off=0, y_off=0, niter=256):
""" A mandelbrot set of geometry (width x height) and iterations 'niter' """
w,h = width, height
pixels = np.arange(w*h, dtype=np.uint16).reshape(h, w)
# The mandelbrot set represents every complex point "c" for which
# the Julia set is connected or every julia set that contains
# the origin (0, 0). Hence we always start with c at the origin.
for x in range(w):
for y in range(h):
# calculate the initial real and imaginary part of z,
# based on the pixel location and zoom and position values
# We use (x-3*w/4) instead of (x-w/2) to fully visualize
# the fractal along the x-axis
zx = 1.5*(x + x_off - 3*w/4)/(0.5*zoom*w)
zy = 1.0*(y + y_off - h/2)/(0.5*zoom*h)
z = complex(zx, zy)
c = complex(0, 0)
for i in range(niter):
if abs(c) > 4: break
# Iterate till the point c is outside
# the circle with radius 2.
# Calculate new positions
c = c**2 + z
color = (i << 21) + (i << 10) + i * 8
pixels[y,x] = color
return pixels
def display(width=1024, height=768, zoom=1.0, x_off=0, y_off=0, cmap='viridis'):
""" Display a mandelbrot set of width `width` and height `height` and zoom `zoom`
and offsets (x_off, y_off) """
pixels = mandelbrot_set(width, height, zoom=zoom, x_off=x_off, y_off=y_off)
# Let us turn off the axes
plt.axis('off')
# to display the created fractal
plt.imshow(pixels, cmap=cmap)
plt.show()
So much for the code. Let us look at some images.
Fractal Images
Here is the beautiful image of the Mandelbrot fractal that our code generates.
Let us play around with the code a bit and try and generate different images.
Here is the fractal set zoomed at 2X with an x offset of -300 pixels.
>>> import mandelbrot
>>> mandelbrot.display(zoom=2.0, x_off=-300)
As you may have expected to see, the fractal is self-similar at higher magnifications.
Here is the fractal with a different color-map than the default one.
>>> import mandelbrot
>>> mandelbrot.display(cmap='magma')
You can also play around with the color scheme to generate differently colored images. For example, by changing the color scheme to $ (i<<4) + (i<<12) + i*2 $, you get the following image, with the green color more accentuated around the edges.
For more information on color scheme experiments with the Mandelbrot Set, visit this discussion.
Concurrent version
It is possible to make this version of the code concurrent and improve its performance by using the multiprocessing
module, just like we did for the Julia set in the 2nd part of this series. It is left as an exercise to the reader.
All the code developed in this post is available in our fractals code repository @ github .
Conclusion
I hope you enjoyed this post in the fractal posts series. In future posts, we will explore the mathematical properties of these fractals and also go on to develop code for more interesting fractals.
Note that name and e-mail are required for posting comments