Homework 4, Problem 1: Creating
the Mandelbrot Set
[40
points; individual or pair]
Submission: Submit your hw4pr1.py file to Canvas
In this lab you will build a
program to visualize and explore the points in and near the Mandelbot Set.
Objectives:
- Use loops and nested loops to solve complex problems (quite
literally!)
- Develop a program using incremental design, i.e., by
starting with a simple task and gradually adding levels of complexity
- Draw connections with mathematics and other disciplines that use
fractal modeling
Introduction to the Mandelbrot Set
The Mandelbrot set is a set of points in the
complex plane that share an interesting property. Choose a complex number c . With this c in mind,
start with
z0 =
0
and then repeatedly iterate as follows:
zn+1 = zn2 + c
The Mandelbrot set is the collection of all complex
numbers c
such that this process does not diverge to infinity as n gets large. There
are other, equivalent definitions of the Mandelbrot set.
The Mandelbrot set is a fractal, meaning that
its boundary is so complex that it cannot be well-approximated by
one-dimensional line segments, regardess of how closely one zooms in on it. In
fact, the Mandelbrot set's boundary has a dimension of 2 (!), though the
details of this are better left to the many available references.
Finding whether a complex number is in
the set or not
Your first task is to write a function named inMSet( c, n ) that takes as input a
complex number c and an integer n. Your
function will return a boolean: True if the complex number c is
in the Mandelbrot set and
False otherwise.
In Python a complex number is represented in terms of
its real part x and its imaginary part y. The
mathematical notation would be x+yi, but in Python the
imaginary unit is typed as 1.0j or 1j, so that
c = x + y*1j
would assign the variable c to
the complex number with real part x and imaginary part y.
Unfortunately, x
+ yj does
not work, because Python thinks you're using a variable named yj.
Also, the value 1
+ j is
not a complex number: Python assumes you mean a variable named j unless
there is an int or a float directly in front of it. Use 1 + 1j instead.
Try it out
Just to get familiar with complex numbers, at the Python prompt
try
>>> c = 3 + 4j
>>> c
(3+4j)
>>> abs(c)
5.0
>>> c**2
(-7+24j)
Python is happy to use the power operator (**)
and other operators with complex numbers. However, note that you cannot compare
complex numbers directly (they do not have an ordering defined on them since
they are essentially points in a 2D space). So you cannot do something like c > 2. However, you CAN compare the magnitude (i.e.,
the absolute value) of a complex number to the number 2.
Back
to the inMSet function...
To determine whether or not a number c is in the
Mandelbrot set, you simply need to start with z0 = 0 + 0j and repeatedly iterate zn+1 = zn2 + c to see if
this sequence of zs stays bounded, i.e., the magnitude of z does not go
to infinity.
Determining whether or not this sequence really goes
off to infinity would take forever! To check this computationally, we will have
to decide on two things:
- the number of
times we are willing to wait for the zn+1 = zn2 + c
process to run
- a value that
will represent "infinity"
The first value above (the number of times we are
willing to run the z-updating process) is the second input, n, to
the function inMSet( c, n ). This is a value you will
want to experiment with, but you might start with 25 as a reasonable initial
value.
The second value above, one that represents infinity,
can be surprisingly low! It has been shown that if the absolute value of the
complex number z
ever gets larger than 2 during its update process, then the
sequence will definitely diverge to infinity. There is no equivalent rule that
tells us that the sequence definitely will not diverge, but it is very likely it
will stay bounded if abs(z)
does not exceed 2 after a reasonable number of iterations, and n is
that "reasonable" number.
So, to get started, write the function inMSet( c, n ) that takes in a complex
number c
and a number of iterations n. It should output False if the sequence zn+1 = zn2 + c ever
yields a z whose magnitude is greater than 2. It returns True otherwise.
Check your inMSet
function with these examples:
>>> n = 25
>>> c = 0 + 0j # this
one is in the set
>>> inMSet( c, n )
True
>>> c = 0.3 + -0.5j #
this one is also in the set
>>> inMSet( c, n )
True
>>> c = -0.7 + 0.3j #
this one is NOT in the set
>>> inMSet( c, n )
False
>>> c = -0.9+0.4j #
also not in the set
>>> inMSet( c, n )
False
>>> c = 0.42 + 0.2j
>>> inMSet( c, 25 ) #
this one seems to be in the set
True
>>> inMSet( c, 50 ) # but
on further investigation, maybe not...
False
Getting
too many Trues?
If so, you might be checking for abs(z) > 2 after the
for loop finishes. There is a subtle reason this won't work. Many values get so
large so fast that they overflow the capacity of Python's floating-point
numbers. When they do, they cease to obey greater-than / less-than
relationships, and so the test will fail. The solution is to check whether
the magnitude of z ever gets bigger than 2 inside the for
loop, in which case you should immediately return False. The return True,
however needs to stay outside the loop!
As the last example
illustrates, when numbers are close to the boundary of the Mandelbrot set, many
additional iterations may be needed to determine whether they escape. We will
not worry too much about this, however.
In order to start creating images with Python,
download
the bmp.py
file at this link
to
the same directory where your hw4pr1.py file is located. You will
want to have the folder holding those files open, because it is that folder (or
Desktop) where the images will be created.
Here is a function to get you started -- try it out!
from bmp import *
def testImage():
""" a function to demonstrate
how
to create and save a bitmap image
"""
width = 200
height = 200
image = BitMap( width, height )
# create a loop in order to draw some pixels
for col in range(width):
if col % 10 == 0: print('col
is', col)
for row in range(height):
if col % 10 == 0 and row %
10 == 0:
image.plotPoint( col, row )
# we have now looped through every image
pixel
# next, we write it out to a file
image.saveFile( "test.bmp" )
Save this function, and then run it by typing testImage(), with the parentheses, at
the Python shell.
If everything goes well, a few print statements will
run, but nothing else will appear. However, if it worked, this function will
have created an image in the same folder where bmp.py and hw4pr1.py are located. It will be called test.bmp and will be a bitmap-type image.
Both Windows and Mac computers have nice built-in
facilities for looking at bitmap-type images. Simply double click on the icon
of the test.bmp image, and you will see
it. For the above function, it should be all white except for a regular, sparse
point field, plotted wherever the row number and column number were both
multiples of 10.
You can zoom in and out of bitmaps with the built-in
buttons on Windows; on Macs, the commands are under the view menu.
Before changing the above code,
write a short comment or triple-quoted string in your hw4pr1.py file describing what would happen if you changed the line
if col % 10 == 0 and row % 10 == 0:
to the line
if col % 10 == 0 or row % 10 == 0:
Then, make that change from and to or and try it.
Just for practice, you might try making larger images -
or images with different patterns - by changing the testImage function appropriately.
Some notes on how
the testImage function works
There are three lines of the testImage function
that warrant a closer look:
- image = BitMap( width, height )
This line of
code creates an image of type BitMap with the
specified height and width. The image variable
holds the whole image! This is similar to the way a single
variable, often called L can hold an arbitrarily
large list of items. When information is gathered together into a list or
image or another structure, it is called a software object or just
an object. We will build objects of our own in a couple of weeks,
so this is a good opportunity to use them without having to worry about
creating them.
- image.plotPoint( col, row ) An important property of software objects
is that they can carry around and call functions of their own! They do
this using the dot . operator. Here, the image object is calling its own plotPoint function in order to, well, plot a point at the
given column and row. Functions called in this way are sometimes called methods.
- image.saveFile( "test.bmp" ) This line actually creates the new test.bmp file that holds the bitmap image. It demonstrates
another method (function) of the software object named image.
From pixel to complex coordinates
Ultimately, what we are trying to do is plot the
Mandelbrot set on a complex coordinate system. However, when we plot points on
the screen, we manipulate pixels. As the example above shows, pixel values
always start at (0, 0) (in the lower left) and grow to (width-1, height-1) in
the upper right. In the example above width and height were both 200, giving us a reasonable size
image. However, if we created a one-to-one mapping between pixel values and
points in the complex plane our real value would be forced to range from 0 to
199 and our imaginary value would be forced to range between 0 and 199. This is
not a very interesting region to plot the Mandelbrot set, because the vast
majority of these points are not in the set (so our image would be mostly
empty).
What we would like to do is to be able to define any
complex coordinate system we want and visualize the points in that coordinate
system. In order to do this, we need to project the pixel values (i.e., the
location were we will possibly put a dot on the screen) onto this coordinate
frame (i.e., the actual
points in the complex plane).
In other words, given a pixel, in order to determine whether to put a dot
there, you need to first figure out what point it corresponds to on the complex
plane and then test to see whether that complex point is in the Mandelbrot set.
To help with this conversion back and forth from pixel
to complex coordinates, write a function named scale( pix, pixNum, floatMin,
floatMax ) that
takes in four things: pix, an
integer representing a pixel column, pixNum, the total number of pixel
columns available, floatMin, the floating-point lower
endpoint of the image's real axis (x-axis), and floatMax, the
floating-point upper endpoint of the image's real axis (x-axis).
The idea is that scale will return the
floating-point value between floatMin and floatMax that corresponds to the position of the
pixel pix, which is somewhere between
0 and pixNum. This diagram illustrates the geometry of
these values:
Once you have written your scale function, here are some test cases to try
to be sure it is working:
>>> scale(100, 200, -2.0, 1.0 ) # halfway from -2 to 1
-0.5
>>> scale(100, 200, -1.5, 1.5 ) # halfway from -1.5 to 1.5
0.0
>>> scale(100, 300, -2.0, 1.0 ) # 1/3 of the way from -2 to 1
-1.0
>>> scale(25, 300, -2.0, 1.0 )
-1.75
>>> scale(299, 300, -2.0, 1.0 ) # your exact value may differ slightly
0.99
Note Although initially stated
as something to find x-coordinate floating-point values, your scale function works equally well for either the
x- or y- dimension. You won't need a separate function for the other axis!
Visualizing the Mandelbrot set in
black and white
This part asks you to put the pieces from the above
sections together into a function named mset( width, height )
that computes the set of points in the Mandelbrot set on the complex plane and
creates a bitmap of them. For the moment, we will use images of width width
and height height
and we will limit the part of the plane to the ranges
-2.0 ≤ x or
real coordinate ≤ +1.0
and
-1.0 ≤ y or
imaginary coordinate ≤ +1.0
which is a 3.0 x 2.0 rectangle.
How to get started?
Start by copying the code from the testImage
function -- those pixel loops using col and row
will be the basis for your mset function, as well. There are at least
five things that will have to change:
- The width
and height
are inputs, not hand-set values.
- For each pixel
col,
you need to compute the real coordinate of that pixel in the
complex plane. Use scale - a reasonable name for the result is
x.
- For each pixel
row,
you need to compute the imagary coordinate of that pixel in the
complex plane. Use scale - a reasonable name for the result is
y.
Even though this will be the imaginary part of a complex number, it
is simply a normal floating point value.
- Using the real
and imaginary parts computed in the prior two steps, create a variable
named c
that holds a complex value with those real and imaginary
parts, respectively. Recall that you'll need to multiply y*1j,
not y*j!
- Your test for
which pixel col
and row
values to plot will involve inMSet, the first function
you wrote.
The variable width and height
inputs provide the opportunity to create both large images (if you're patient)
and small images (for quicker prototyping to see if things are working). You
might first try
>>> mset( 300, 200 )
to be sure that the image you get is a black-and-white
version of the many Mandelbrot set images on the web. Then, you might want to
try a larger one.
Adding features to your mset function
Magic Numbers
are simply literal numeric values that you have typed into your code. They're
called magic numbers because if someone tries to read your code, the values and
purpose of those numbers seem to have been pulled out of a hat... . For
example, your mset function may have called inMSet( c, 25 ). A newcomer to your
code (and this problem) would have no idea what that 25 represented or why you chose that
value.
To keep your code as flexible and expandible as
possible, it's a good idea to avoid using these "magic numbers" for
important quantities in various places in your functions. Instead, it's a good
idea to collect all of those magic numbers at the very top of your functions
(after the docstring) and to give them useful names that suggest their purpose.
For example, you might have the line
numIterations = 25 #
after this many tries, we assume z in in the MSet
and then the function call would look like
inMSet( c, numIterations )
In addition to being clearer, this makes it much
easier to add functionality to your code -- all of the important quantities are
defined one time and in one place.
For this part of the lab, move
all of your magic numbers to the top of the function and give them descriptive
names. This will help for the next two parts, as well.
Changing the number of iterations
Next, experiment with different values of numIterations -- or whatever you are calling the maximum
number of update steps of z before your inMSet returns True. In particular, you might
try the values
numIterations = 5
numIterations = 10
and
numIterations = 50
in addition to our default value of 25. To make these tests more efficient, you
could make numIterations an input to mset but it's certainly not required.
Zooming
by changing the range
Another advantage of having magic numbers at the top of
a function is that it becomes very easy to make them inputs of the function
instead of static, hand-assigned values.
For this part, change the inputs to mset so that the function takes in a list:
mset( width, height, coordinateList )
where coordinateList
is
a list of four endpoints of the complex plane's axes. For example,
mset( width, height, [ -2.0, 1.0, -1.0, 1.0 ] )
would specify a horizontal (real or x-axis) range from -2.0 to 1.0 and would specify a
vertical (imaginary or y-axis) range from -1.0 to 1.0.
Once you have this coordinateList input working, use it to "zoom in" on an interesting
piece of the Mandelbrot set's boundary.
One difficulty here is that the ratio of height/width
needs to be the same as the (vertical range)/(horizontal range). To get around
this, remove the height input and instead compute
it based on the fact that
height = width*(vertical range)/(horizontal range)
This will ensure that the aspect ratio of the
resulting image is true.
Remember that if you like a particular image and don't
want to overwrite it, you can simply change its name from test.bmp to something else.
Completely Optional Extensions
The next two extensions are completely optional. You
don't have to do them, but they are fun... and the second is worth a little
extra credit.
You don't have to use black and white! The image
object representing a BitMap? supports changes to the
"pen color" through the method call
image.setPenColor( Color.GREEN )
where the actual color can be any of BLACK, BLUE,
BROWN, CYAN, DKBLUE, DKGREEN, DKRED, GRAY, GREEN, MAGENTA, PURPLE, RED, TEAL,
WHITE, or YELLOW . By setting both the pixels out of the set as well as
those in the set, you can create many variations in background and
foreground... .
Try it out!
Visualizing escape velocities in
greyscale or color [+ 5 e.c. points]
Images of fractals often use color to represent how
fast points are diverging toward infinity when they are not contained in
the fractal itself. For this problem, create a new version of mset
called msetColor .
Simply copy-and-paste your old code, because its basic bahavior will be the
same. However, you should alter the msetColor function so that it plots
points that escape to infinity more quickly in different colors.
Thus, have your msetColor function plot the
Mandelbrot set as before. In addition, however, use at least three different
colors to show how quickly points outside the set are diverging -- this is
their "escape velocity." Making this change will require a change to
your inMSet
helper function, as well. Perhaps copy, paste, and rename inMSet
so that you can change the new version without affecting the old one.
There are several ways to measure the "escape
velocity" of a particular point. One is to look at the resulting magnitude
after the iterative updates. Another is to count the number of iterations
required before it "escapes" to a magnitude that is greater then 2.
An advantage of the latter approach is that there are fewer different escape
velocities to deal with.
Choose one of those approaches
-- or design another of your own -- in order to implement msetColor .
Note on defining your own colors
You are not limited to using the predefined colors listed above. You can create
a color object, similar in spirit to the bitmap objects that we have
been working with, through the following code:
mycolor = Color(red,green,blue)
where red, green, and blue
are integers between 0
and 255
inclusive. They indicate how much of those color components the object mycolor
will have. After this definition, mycolor may be used just like any of the
predefined colors.
If you have gotten to this point, you have completed Lab 1! You should submit your hw4pr1.py file at Canvas .
Next