Chapter 2 The Mandelbrot Set
But better die than live mechanically
a life that is a repetition of repetitions.
— D. H. Lawrence, Women in Love
The Mandelbrot set is one of the most fascinating objects in mathematics. It is based on a computational process that repeatedly applies the mapping \(z \mapsto z^2+c\), where \(z\) and \(c\) are complex numbers and the starting value is \(z_0=0\). This process gives us the following infinite sequence \((z_k)\) of complex numbers:
\(\seteqnumber{0}{2.}{0}\)\begin{align*} z_0 &= 0\\ z_1 &= z_0^2 + c\\ z_2 &= z_1^2 + c\\ \dots &= \dots \end{align*} The behavior of this sequence depends on the parameter \(c\). For some values, the sequence diverges, which means that the numbers \(z_k\) gradually move away from the origin. But for other values of \(c\), the sequence is bounded: all numbers \(z_k\) stay in a fixed circle around the origin. If we now examine all possible values of \(c\) and mark those for which the sequence \((z_k)\) is bounded, we obtain the pear-shaped object at the center of Fig. 2.1. This is the Mandelbrot set.
The real beauty of this set emerges if we magnify any region close to the boundary of this black shape and add a little color; four such regions are shown in the figure. No matter where you look or how far you zoom into the boundary, new variations of spirals, tendrils, and ornaments keep appearing. In this chapter we will discuss how to compute these images ourselves.
2.1 Computing the Mandelbrot Set¶
Let’s start with a few preliminaries before we turn to the problem of computing and visualizing the Mandelbrot set. The mapping \(z\mapsto z^2+c\), and therefore the Mandelbrot set itself, is defined in terms of complex numbers. You can think of a complex number as a pair of real numbers that is written in the form \(z=x + \I y\), where \(x\) is called the real part of \(z\) and \(y\) its imaginary part. The symbol \(\I \) is called the imaginary unit. There is a one-to-one correspondence between points \((x,y)\) in the plane and complex numbers \(x+\I {}y\). Complex numbers are therefore often said to live in the complex plane, and this geometric interpretation helps us visualize many of their properties. But compared to geometric points, complex numbers have a much richer algebraic structure: We can add, subtract, multiply, and divide complex numbers, they follow the same laws of algebra as the real numbers, and it is even possible to generalize common functions such as powers, square roots, \(\sin \), and \(\exp \) to complex arguments.
To compute the Mandelbrot set, we will only need a few elementary arithmetic operations. To add and subtract complex numbers, you simply add or subtract their real and imaginary parts separately. For two complex numbers \(z = x+\I y\) and \(w = u+\I v\) we have:
\(\seteqnumber{0}{2.}{0}\)\begin{align} z + w &= (x + u) + \I (y + v)\\ z - w &= (x - u) + \I (y - v) \end{align} The rules for multiplication and division are slightly more complicated and less obvious. The product of two complex numbers is defined as
\(\seteqnumber{0}{2.}{2}\)\begin{equation} z w = (x u - y v) + \I (x v + y u) \end{equation}
and their quotient as
\(\seteqnumber{0}{2.}{3}\)\begin{equation} z / w = \frac {x u + y v}{u^2 + v^2} + \I \frac {y u - x v}{u^2 + v^2}. \end{equation}
These definitions ensure that most of the familiar laws from real arithmetic carry over to complex numbers and that expressions involving complex numbers can be manipulated as normal. For instance, complex addition and multiplication are commutative and associative, the distributive law \((a+b)c = ac + bc\) holds, and we have \(a/b = a\cdot 1/b\).
The most famous property of complex numbers is the fact that the square of the imaginary unit is \(-1\), which follows directly from the definition of multiplication:
\(\seteqnumber{0}{2.}{4}\)\begin{equation*} \I ^2 = (0+1\I )(0+1\I ) = -1 + 0\I = -1. \end{equation*}
Thanks to this identity, the square root function can be extended to negative numbers when working with complex numbers. For instance, we have \(\sqrt {-1}=\I \) and \(\sqrt {-4}=2\sqrt {-1}=2\I \). The absolute value (or modulus) of a complex number \(z=x+\I y\) is defined as its distance from the origin: \(|z|=\sqrt {x^2 + y^2}\).
The Mandelbrot Iteration
The Mandelbrot set is based on the following infinite sequence of complex numbers \(z_k\) that starts at the origin:
\(\seteqnumber{0}{2.}{4}\)\begin{equation} \label {eq:mandelbrot} \begin{aligned} z_0 &= 0,\\ z_k &= z_{k-1}^2 + c. \end {aligned} \end{equation}
The parameter \(c\) is an arbitrary complex number, so these equations define an endless number of infinite sequences. We also write \(z_k(c)\) to emphasize the dependence of the sequence on the parameter \(c\).
Let’s take a look at a few examples to see how these sequences behave. The simplest case is \(c=0\), which gives us \(z_k=0\) for every \(k\); the resulting sequence consists of an infinite number of zeros \(0,0,0,\dots \). For \(c=1\) we obtain the sequence \(0, 1, 2, 5, 26, 677, \dots \) in which each number is \(1\) plus the square of the previous number; this sequence grows without bound. Finally, for the complex unit \(c=\I \), we obtain a sequence that starts with the numbers \(0\), \(\I \), \(-1 + \I \), \(-\I \), \(-1+\I \) and then continue to alternate between the last two numbers.
As mentioned in the introduction, the Mandelbrot set consists of those numbers \(c\) for which the sequence \(z_k(c)\) stays bounded, that is, for which all numbers in the sequence \(z_k\) remain inside some finite region around the origin. We just saw that the sequences for \(c=0\) and \(c=\I \) are bounded, so these two numbers are part of the Mandelbrot set, but the number \(c=1\) is not in the set because its sequence diverges.
We can now draw a black-and-white image of the Mandelbrot set by choosing a region of the complex plane, dividing it into width times height square “pixels” and setting each one to either black or white depending on whether Mandelbrot iteration for the complex number at its center produces a bounded or an unbounded sequence (Algorithm 2.1).
But how do we determine whether, for a given value of \(c\), the sequence \(z_k(c)\) is bounded or unbounded? The goal is to develop a criterion that is based on computing and analyzing a finite number of terms \(z_0, z_1, \dots , z_n\) and that gives us the correct answer with a reasonable level of accuracy.
Fortunately, the Mandelbrot iteration has an important property that makes it easy to identify the majority of unbounded sequences. The key insight is that the sequence \(z_k(c)\) is guaranteed to diverge if there is just one \(j\) for which \(|z_j(c)|>2\). Conversely, every bounded sequence must satisfy \(|z_k(c)|<2\) for all \(k\), that is, the sequence in its entirety is confined to a circle of radius 2. The number 2 is also called the escape radius of the Mandelbrot iteration because the sequence diverges or “escapes” once it crosses this threshold. A proof of this result is provided in the following aside.
Box 2.1: Escape radius of the Mandelbrot iteration Here is the proof that the Mandelbrot sequence \(z_k(c)\) diverges once it crosses the escape radius 2. We first show that under certain conditions every step of the Mandelbrot iteration moves us further away from the origin, and we then prove that these conditions are always met.
First, assume that the \(j\)th element of the sequence \(w=z_j(c)\) satisfies both \(|w|>2\) and \(|w|>|c|\). If this is the case, every additional step of the Mandelbrot iteration takes us even further away from the origin because the absolute value \(|w|\) is less than the absolute value of the next element in the sequence \(|w^2+c|\):
\(\seteqnumber{0}{2.}{5}\)\begin{align*} \frac {|w^2+c|}{|w|} &\ge \frac {|w|^2 - |c|}{|w|} && \text {triangle inequality}\\ &= |w| - \frac {|c|}{|w|}\\ &> |w| - 1 && \text {since $|w|>|c|$}\\ &> 1 && \text {since $|w|>2$} \end{align*}
We now show that every sequence that crosses the escape radius must have an element \(w\) that satisfies both \(|w|>2\) and \(|w|>c\). We distinguish two cases. If \(|c|\le 2\), the first \(z_m\) that lies outside the escape radius obviously satisfies \(|z_m|>2\ge |c|\), so we are already done. Otherwise, if \(|c|>2\), we consider the first three elements of the sequence:
\(\seteqnumber{0}{2.}{5}\)\begin{equation*} z_0=0,\quad z_1 = c,\quad z_2=c^2+c. \end{equation*}
It’s easy to see that the element \(z_2\) satisfies the conditions above since
\(\seteqnumber{0}{2.}{5}\)\begin{equation*} |z_2| = |c^2+c|\ge |c^2|-|c| = |c|(|c|-1) > |c|\ge 2. \end{equation*}
Therefore, regardless of the value of \(c\), the Mandelbrot iteration diverges once it crosses the escape radius.
Let us define the escape time of a number \(c\) as the smallest index \(j\) for which \(z_j(c)\) lies outside the escape radius. If the sequence is bounded, it must stay inside the escape radius indefinitely, so its escape time is infinite. The relationship beteween the escape radius and the escape time is illustrated in Fig. 2.2, which shows the Mandelbrot set embedded inside a circle of radius 2. The three paths shown in the image illustrate the behavior of sequences \(z_k(c)\) for different parameters \(c\). Two of the paths remain inside the escape radius indefinitely, so their escape time is infinite and the respective starting points belong to the Mandelbrot set. The third path, on the other hand, leaves the circle after five steps and then escapes to infinity. The escape time of the starting point \(c=-0.9+0.5\I \) is therefore \(5\) and the number does not belong to the Mandelbrot set.
To compute the escape time of a number \(c\), we first compute successive values \(z_0, z_1, \dots \) of the Mandelbrot iteration. If we reach a number with \(|z_j|>2\) we stop immediately and return \(j\) as the escape time. Otherwise, we continue until we reach a certain iteration limit MaxIter. Every sequence that is still inside the escape radius after MaxIter steps is considered bounded. The larger we choose the iteration limit, the more accurate our computation becomes. A small iteration limit of \(100\) is sufficient for drawing low-resolution images of the Mandelbrot set, but larger values are necessary to resolve fine details at higher magnifications.
The escape time is also the key to creating images of the Mandelbrot set that are far more rich and detailed than the simple black-and-white images we used so far. Figure 2.3 shows the result of mapping different escape times to different colors. Numbers outside the escape radius (escape time 0) are shown as light blue pixels, and as we get closer to the boundary of the Mandelbrot set, the escape time increases and slowly shifts from blue to red and yellow (escape time 80). Numbers with an infinite escape time are still shown as black pixels. The algorithm for drawing the Mandelbrot set in color is summarized in Algorithm 2.2.
Exercises
[2.1]The best way to appreciate the beauty and complexity of the Mandelbrot set is to explore it interactively. Use one of the many free fractal viewers available on the Internet to find your favorite spots in the Mandelbrot set. Here are the coordinates of a few interesting regions to get you started: the seahorse valley around \(-0.75+0.1\I \); the elephant valley around \(0.3 + 0\I \); and the northwest bulb around \(-0.5434 + 0.6132\I \).
2.2 Drawing the Mandelbrot Set¶
Now that we understand the basic algorithm for computing the Mandelbrot set, we can turn to its implementation. Let’s start with the Mandelbrot iteration
\(\seteqnumber{0}{2.}{5}\)\begin{equation*} z_0 = 0,\qquad z_k = z_{k-1}^2+c, \end{equation*}
which computes a sequence of complex numbers. Instead of defining a custom data type for representing complex numbers we take the slightly simpler route and implement the necessary operations directly. To do this, we first split the numbers \(z_k\) and \(c\) into their real and imaginary parts by writing \(z_k=x_k + \I y_k\) and \(c=c_x+\I c_y\) and then expand the expression \(z_{k-1}^2+c\) using the rules of complex arithmetic. This gives us two new sequences, one for the real parts \(x_0,x_1,\ldots \) and one for the imaginary parts \(y_0,y_1,\ldots \):
\(\seteqnumber{0}{2.}{5}\)\begin{equation} \label {eq:mandelbrot:real} \begin{aligned} x_0 &= 0 & x_k &= x_{k-1}^2 - y_{k-1}^2 + c_x && \text {real part},\\ y_0 &= 0 & y_k &= 2x_{k-1} y_{k-1} + c_y && \text {imaginary part}. \end {aligned} \end{equation}
Since all numbers and variables that occur in Eq. (2.6) are real-valued, the expressions can be translated directly into Java code that operates on floating-point numbers.1
We start in Listing 2.1 with the implementation of a function escapeTime() that computes the escape time of a complex number \(c\). The first two arguments cx and cy are the real and imaginary parts of \(c\), and the third argument maxIter is the iteration limit. The return value of escapeTime() is either the number of iterations at which the trajectory leaves the escape radius or maxIter if the trajectory doesn’t escape. To check whether the complex number \(z_k\) lies outside the escape radius, the function compares the squared length \(|z_k|^2=x_k^2+y_k^2\) to the square of the escape radius. Comparing squared lengths (or distances) instead of regular lengths is a common programming trick that saves one square root per comparison. Especially for programs that would otherwise have to compute a large number of square roots, this is a simple but effective optimization.
Listing 2.1ch2β―/β―MandelbrotSet
public static int escapeTime(double cx, double cy, int maxIter) {
double escapeRadius = 2.0;
double xn = 0.0, yn = 0.0;
for (int k = 0; k < maxIter; k++) {
if (xn * xn + yn * yn > escapeRadius * escapeRadius)
return k;
double x = xn, y = yn;
xn = x * x - y * y + cx;
yn = 2 * x * y + cy;
}
return maxIter;
}
We can now draw an image of the Mandelbrot set by iterating over each pixel and setting its color based on the escape time of the corresponding complex number. Suppose we want to draw the Mandelbrot set in a certain rectangular region of the complex plane and produce a raster image that is \(w\) pixels wide and \(h\) pixels high. Figure 2.4 illustrates how we translate between pixels and complex numbers: We first divide the region into a grid of \(w\times h\) pixels, and then compute the coordinates of the complex number at the center of each square.
The drawMandelbrot() function shown in Listing 2.2 implements these computations. The function takes three arguments: The first is the region of the complex plane to be drawn, the second argument colors[] specifies the color palette that maps escape times to colors, and the last is the image we are drawing into.2 The function starts by computing the size of each pixel and the coordinates xMin and yMin of the complex number at the center of the lower-left pixel. We assume that the array colors specifies one color for each possible escape time, so we set the iteration limit maxIter to the index of the last entry in colors; this last entry holds the color that is used for the interior of the Mandelbrot set and is conventionally set to black. The function then iterates over the pixels of image and colors them according to the escape time of the corresponding complex number.
Listing 2.2ch2β―/β―MandelbrotPainter
public static void drawMandelbrot(Rectangle region, int[] colors,
BufferedImage image) {
double pixelWidth = region.width() / image.getWidth();
double pixelHeight = region.height() / image.getHeight();
double xMin = region.bottomLeft().x() + pixelWidth / 2;
double yMin = region.bottomLeft().y() + pixelHeight / 2;
int maxIter = colors.length - 1;
for (int y = 0; y < image.getHeight(); y++) {
for (int x = 0; x < image.getWidth(); x++) {
int escapeTime = MandelbrotSet.escapeTime(
xMin + x * pixelWidth, yMin + y * pixelHeight,
maxIter);
image.setRGB(x, image.getHeight() - 1 - y,
colors[escapeTime]);
}
}
}
JDK: BufferedImage
•MandelbrotSet: escapeTime()
•Rectangle: bottomLeft(), height(), width()
Finally, let us use drawMandelbrot() to create a series of images that zooms into a particular region of the Mandelbrot set, as shown in Fig. 2.5. Each image magnifies the central region of the previous image by a factor of 2.5, so the total magnification of the last image in the lower-right corner is \(2.5^{14}\approx 370\,000\).
To create such a sequence of images, we call drawMandelbrot() multiple times and in each step shrink the region that is being drawn by a constant factor while keeping the center fixed. The corresponding rectangles in the complex plane can be computed as shown in Listing 2.3. The function first computes the width and height of the rectangle in the complex domain. The “magic constant” 2.25 used in the computation of width is chosen so that a magnification of 1 produces a rectangle that encloses the entire width of the Mandelbrot set. The height of the of the region is then obtained by scaling width by the aspect ratio of the resulting image, that is, the ratio of its width to its height. We can then easily compute the lower-left and upper-right corners of the region in the complex plane, which are used to initialize the Rectangle returned from the function.
Listing 2.3ch2β―/β―MandelbrotPainter
public static Rectangle zoomRect(Point center,
double magnification,
int imageWidth, int imageHeight) {
double w = 2.25 / magnification;
double h = w / imageWidth * imageHeight;
return new Rectangle(
new Point(center.x() - w / 2, center.y() - h / 2),
new Point(center.x() + w / 2, center.y() + h / 2));
}
We can now create a series of images that keep the center fixed and repeatedly increase the magnification by a constant factor:
var center = new Point(-0.5364866407792761, 0.5230664066236846);
double magnification = 1.0;
for (int i = 0; i < 10; i++) {
var region = zoomRect(center, magnification,
image.getWidth(), image.getHeight());
drawMandelbrot(region, colors, image);
// Display or save image
magnification *= 2.5;
}
•MandelbrotPainter: drawMandelbrot(), zoomRect()
•Point: Point()
1 Exercise 2.2 asks you to derive these equations yourself.
2 The Rectangle type was introduced in Exercise 1.4.
Exercises
[2.2]Use the rules of complex arithmetic to show that the two versions of the Mandelbrot iteration in (2.5) and (2.6) are equivalent.
[2.3]The Mandelbrot set can be rendered in countless different ways by changing the mapping from escape times to colors. Create images with different color mappings (\(k\) is the escape time):
-
1. Black and white: black if the number belongs to the Mandelbrot set and white otherwise.
-
2. Grayscale: either increase the intensity linearly with \(k\) or use an oscillating function such as \(0.5 + 0.5\sin (2k)\).
-
3. Color: for each color channel, vary the intensity using the function \(0.6 + 0.4\cdot \cos (f_c\cdot k + \phi _c)\), where \(f_c\) and \(\phi _c\) denote independent parameters for the red, green, and blue channels. (The colored images in this chapter were generated using \(f_{rgb}=(0.05, 0.04, 0.03)\), and \(\phi _{rgb} = (2.8, 1.7, 0.5)\).)
2.3 The Buddhabrot¶
As a final example of iteration and fractals, let’s take a look at a particularly interesting method for visualizing the Mandelbrot set that visualizes not just the escape time but the way the Mandelbrot iteration distributes points over the complex plane. The result is the Buddhabrot, shown in Fig. 2.6, which owes its name to the fact that the resulting shape resembles a sitting Buddha.3
To draw the Buddhabrot, we focus our attention on unbounded trajectories of the Mandelbrot iteration in Eq. (2.5). Recall that unbounded trajectories are characterized by the fact that they eventually leave the escape radius, but before they disappear from view, they visit a finite set of numbers inside the escape radius. The Buddhabrot shows which numbers are, on average, visited most frequently: The complex numbers in the bright regions of Fig. 2.6 are traversed more frequently than those in the dark regions.
To compute the Buddhabrot, we divide the region of the complex plane we want to display into pixels as before, but this time we associate a counter with each pixel. We then generate a random complex number inside the escape radius and compute its the trajectory using the normal Mandelbrot iteration. If the trajectory is bounded we ignore it, but if it is unbounded we take each complex number on the trajectory, determine the corresponding pixel, and increase its counter. If we repeat this process for many random trajectories, we obtain a histogram that reflects how complex numbers are distributed by the Mandelbrot iteration.
Since the code for computing the Buddhabrot is slightly more complex than the one for the regular Mandelbrot set, we collect all the required data and functions in a class called Buddhabrot (Listing 2.4). As before, the imageWidth and imageHeight fields hold the width and height of resulting image, and region is the region in the complex plane that is being displayed. The constructor Buddhabrot() initializes imageWidth and imageHeight directly and computes a suitable region of interest using the zoomRect() method, which we already used to compute a region of interest for the Mandelbrot set.
Listing 2.4ch2β―/β―Buddhabrot.java
// A class for computing the Buddhabrot.
public class Buddhabrot {
private final int imageWidth; // width of final image
private final int imageHeight; // height of final image
private final Rectangle region; // region of interest
private Buddhabrot(int imageWidth, int imageHeight) {
this.imageWidth = imageWidth;
this.imageHeight = imageHeight;
region = MandelbrotPainter.zoomRect(
new Point(-0.5, 0.0), 1.0,
imageWidth, imageHeight);
}
// ...
}
•MandelbrotPainter: zoomRect()
•Rectangle
To compute an image of the Buddhabrot, we generate millions of random trajectories and count how often they traverse each pixel. To store the result we use a two-dimensional array of integers that is as large as the final image:
int[][] histogram = new int[imageHeight][imageWidth];
The followTrajectory() function defined in the Listing 2.5 computes a single trajectory that starts at the complex number \(c\) and increment the counters along its path. The arguments cx and cy are the real and imaginary part of \(c\), and maxIter is the number of iterations we want to perform. The function starts by computing the size of a single pixel pixelSize and then performs the Mandelbrot iteration to compute the complex numbers xn, yn along the trajectory. For every such number that lies inside the region of interest, the corresponding histogram entry is incremented; numbers outside of region are ignored.
Listing 2.5ch2β―/β―Buddhabrot
// Increment histogram entries for trajectory starting at cx, cy.
private void followTrajectory(double cx, double cy, int maxIter,
int[][] histogram) {
double pixelSize = region.width() / imageWidth;
double xn = 0.0, yn = 0.0;
for (int k = 0; k < maxIter; k++) {
double x = xn, y = yn;
xn = x * x - y * y + cx;
yn = 2 * x * y + cy;
int row = (int) ((yn - region.bottomLeft().y())
/ pixelSize);
int col = (int) ((xn - region.bottomLeft().x())
/ pixelSize);
if (col >= 0 && col < imageWidth
&& row >= 0 && row < imageHeight)
histogram[row][col]++;
}
}
Each trajectory should start at a random point inside the escape radius. This presents us with the next problem: Most random number generators only produce real numbers in an interval, but we need to generate random points inside a circle. How can this be done?
Let’s first discuss the solution for the unit circle, the circle of radius 1 that is centered at the origin. Every point inside this circle can be specified using polar coordinates \((r,\phi )\), where \(r\) is the distance of the point from the origin and \(\phi \) the angle it makes with the x-axis.4 A random point can then be generated by creating a random number \(u\) between 0 and 1 and a random number \(\phi \) between 0 and \(2\pi \), and then converting back to Cartesian coordinates
\(\seteqnumber{0}{2.}{6}\)\begin{equation} x = \sqrt {u}\cos (\phi ),\qquad y = \sqrt {u}\sin (\phi ). \end{equation}
The reason why the radius must be \(\sqrt {u}\) instead of simply \(u\) is easy to explain: The area of a circular disk grows quadratically with its radius, so the likelihood that a random point occurs at a certain distance from the origin must also increase quadratically with that distance. When generating a random radius, we must therefore transform the uniform value \(u\) using the inverse function, the square root. Random points in circles of arbitrary size can be obtained by multiplying \(x\) and \(y\) by the desired radius.
The function shown in Listing 2.6 uses this technique to create a total of numTrajectories trajectories that start at random points inside the escape radius. Since the Buddhabrot uses only trajectories that escape after a finite number of steps, we update the histogram using followTrajectory() only if the escape time is less than maxIter.
Listing 2.6ch2β―/β―Buddhabrot
// Compute histogram of n random trajectories.
int[][] computeHistogram(int numTrajectories, int maxIter) {
var histogram = new int[imageHeight][imageWidth];
double escapeRadius = 2.0;
for (int i = 0; i < numTrajectories; i++) {
double angle = 2 * Math.PI * Math.random();
double radius = escapeRadius * Math.sqrt(Math.random());
double cx = radius * Math.cos(angle);
double cy = radius * Math.sin(angle);
int escapeTime = MandelbrotSet.escapeTime(cx, cy,
maxIter);
if (escapeTime < maxIter) // unbounded?
followTrajectory(cx, cy, escapeTime, histogram);
}
return histogram;
}
JDK: Math
•Buddhabrot: followTrajectory()
•MandelbrotSet: escapeTime()
Now that we know how to compute the histogram of trajectories, we can draw an image of the Buddhabrot by mapping each histogram entry to a suitable color (Listing 2.7). Similar to the way we colorized the Mandelbrot set, we map each histogram entry to a color from a fixed array of colors. A small but important difference is that the numerical range of the histogram entries is not known in advance, so we first determine the largest entry in the histogram maxValue and use it to reduce every entry to a number in the range \([0,1]\).
Listing 2.7ch2β―/β―Buddhabrot
public static void drawBuddhabrot(int numTrajectories,
int maxIter,
BufferedImage image, int[] colors) {
var buddhabrot = new Buddhabrot(image.getWidth(),
image.getHeight());
var histogram = buddhabrot.computeHistogram(numTrajectories,
maxIter);
int maxValue = maximum(histogram);
for (int y = 0; y < image.getHeight(); y++) {
for (int x = 0; x < image.getWidth(); x++) {
double v = (double) histogram[y][x] / maxValue;
image.setRGB(x, image.getHeight() - 1 - y,
colors[(int) (v * (colors.length - 1))]);
}
}
}
static int maximum(int[][] array) {
int max = Integer.MIN_VALUE;
for (int[] row : array) {
for (int value : row)
max = Math.max(max, value);
}
return max;
}
3 The Buddhabrot visualization is originally due to Melinda Green (http://superliminal.com/fractals/bbrot/bbrot.htm).
4 We discussed polar coordinates in Exercise 1.6.
Exercises
[2.4]The Nebulabrot is an alternative rendering of the Buddhabrot that resembles an intergalactic nebula. The image is obtained by writing three versions of the Buddhabrot to separate color channels of the final RGB image: For the red channel, the iteration limit is set to 5000, for the green channel to 500, and for the blue channel to 50. Write a program that computes the Nebulabrot.
2.4 Problems¶
Drawing and Animating Julia Sets
[2.5]Julia sets are a family of fractals that are closely related to the Mandelbrot set. Both the Mandelbrot set and Julia sets are defined by repeated application of \(z\mapsto z^2+c\), but they differ in which variable is kept constant: For the Mandelbrot set, we use a fixed starting point \(z_0=0\) and study the escape time as a a function of the parameter \(c\), whereas for Julia sets we keep \(c\) constant and study the escape time as a function of the starting value \(z_0\). In contrast to the Mandelbrot set, which is unique, there are infinitely many Julia sets, one for every possible choice of \(c\); two examples are shown in Fig. 2.7.
Julia sets can be visualized by dividing a region of the complex plane into pixels, taking the center of each pixel as the starting point \(z_0\), and coloring the pixel according to the escape time of the trajectory starting at \(z_0\). Again, numbers with an infinite escape time are said to lie inside the Julia set and are usually colored black, and all other pixels are colored by mapping their escape time to a color palette.
-
a) Write a program that computes and displays Julia sets.
Since Julia sets depend on a continuous parameter \(z_0\), they can be animated by moving \(z_0\) along a line, a circle, or another smooth curve and creating a movie from the resulting images.
-
b) Extend your program to generate animations of the Julia set, for example by generating a sequence of images.
The Game of Life
[2.6] Fractals are not the only example of simple rules that give rise to complex behavior. Another example is the Game of Life, or simply Life. The game is played on a rectangular grid of cells, each of which represents a simple life form that can be either dead or alive. The Game of Life starts with an initial population of cells and simulates its evolution over multiple generations. Each cell is influenced by its eight direct neighbors
The evolution of cells over time is governed by the following rules:
-
• A living cell survives if there are exactly two or three other living cells in its neighborhood.
-
• A living cell dies from isolation if there fewer than two living cells in its neighborhood and from overpopulation if there are more than four living cells.
-
• A dead cell becomes alive if there are exactly three living cells in its neighborhood.
Even though the evolution of individual cells in Life is governed by a few simple rules, the behavior of multiple interacting cells can be surprisingly complex. We call an arrangement of cells with interesting properties a life form. One of the simplest life forms is the blinker, which consists of three active cells in a row, and it continually alternates between a horizontal and a vertical configuration:
Another important life form is the glider, which consists of five cells that slowly move across the board:
Every four iterations the glider advances one step to the right and one step down. (Gliders that move along the three other diagonals can be obtained by rotating the initial configuration of cells.) Life forms that move on their own are called spaceships, and the glider is the smallest possible spaceship in Life.
In theory, the grid of cells extends in all directions and life forms can grow infinitely large, but when implementing the Game of Life on a computer, we need to store the game’s state in a finite amount of memory. One common approach is to simulate an infinite board by periodically repeating a finite board in all directions. With such a periodic board, a glider that crosses one edge reappears at the opposite edge.
-
a) Implement the Game of Life. Choose a data structure for representing the grid of cells and implement methods for reading and writing the state of individual cells. Make the board periodic.
-
c) Consider all the life forms that consist of \(n\) living cells in a row. For \(n=3\) you obtain the blinker we discussed above. Study and describe the behavior of the other life forms with \(n\) between 1 and \(10\).
One of the most important discoveries made about the Game of Life is that it is Turing-complete, which means we can construct life forms that can perform any kind of computation — at least theoretically. The derivation of this fascinating result is non-trivial, but it’s easy to see how simple logical functions can be implemented in Life. One of the main ideas is to use streams of gliders to encode and exchange information between different parts of the “computer” and use specialized life forms to operate on these gliders.
-
d) Figure 2.8 shows an arrangement of four gliders that are advancing diagonally from the top-right corner, and a device called a glider gun that is placed in the lower-left corner. Study the evolution of the configuration of cells shown in Fig. 2.8 over the first few hundred generations in a sufficiently large domain. What does the Game of Life “compute” for this input?
2.5 Chapter Notes¶
The first fractals were discovered at the end of the 19th century, and their unexpected and puzzling properties confounded contemporary mathematicians. One such fractal is the Weierstrass function shown in Fig. 2.9, which was discovered in 1872 by the German mathematician Karl Weierstrass. Regardless of where you look or how much you magnify the function’s graph, the curve is always contiguous but it isn’t smooth and it is impossible to compute the function’s slope at any point — in mathematical terms, it is everywhere continuous but nowhere differentiable.
The discovery of the Weierstrass function wasn’t particularly well received: “I turn away with fright and horror of this terrible scourge of continuous functions without derivative,” the British mathematician Charles Hermite wrote of Weierstrass’ function in 1893. However, as more examples of such “terrible scourges” were discovered over the following decades, mathematicians slowly began to appreciate their unusual behavior.
The topic attracted considerable attention in the second half of the 20th century, when electronic computers became powerful enough to perform the laborious computations that underlie all fractals and to visualize the resulting shapes. We now know countless computational processes that produce fractal structures, and they all involve some form of repetition [18]. The term fractal was coined by Benoit Mandelbrot in 1975. It was also Mandelbrot who recognized that fractals — far from being the mathematical aberrations they were considered previously — naturally occur in the study of a wide range of complex systems, from biology and geography to astronomy and economics [68].
For an introduction to the Mandelbrot set and some of its properties, see the article by Dewdney [29]. The Mandelbrot set and Julia sets are also discussed in most books on fractals including the one by Peitgen et al. [76]. An overview of different algorithms for visualizing Mandelbrot sets can be found at https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set.
The process of converting an abstract mathematical description of a scene into an array of colored pixels is called rasterization. In our code for drawing the Mandelbrot set, we computed each pixel’s color from the escape time of a single complex number at its center. If you think of a pixel as a “little square,” you may be wondering whether this is the correct approach: Shouldn’t a pixel’s color be influenced by all complex numbers in its area? The short answer is: it depends [16]. Taking a single sample per pixel is indeed sufficient in regions where the escape time varies slowly. But in regions where it varies quickly, especially close to the boundary of the Mandelbrot set, visual artifacts can occur that usually manifest themselves as jagged boundaries and uneven color transitions. These artifacts are a common problem in image processing and computer graphics known as aliasing. For a discussion of sampling and aliasing in the context of computer graphics, see the book by Pharr et al. [77, Chapter 7]; even though the book’s focus in on ray tracing, the basic ideas apply to image generation in general.
The Game of Life discussed in Exercise 2.6 was invented by John Horton H. Conway in 1970 and quickly became popular in the following years, in large part due to Martin Gardner’s article in Scientific American [42]. Over the years, thousands of different life forms have been discovered and studied in detail; for an overview, see the LifeWiki at https://www.conwaylife.com/wiki/Main_Page. One of the most comprehensive accounts of Life is the book by Johnston and Green [55], which also includes a proof that Life is Turing-complete.
Life belongs to a more general class of computational systems known as cellular automata. A cellular automaton consists of a grid of cells that can change their state in response to other cells in their neighborhood. The simplest cellular automata are rows of binary cells that interact with their left and right neighbors; we will discuss these elementary cellular automata in Exercise 5.15. As in the case of Life, cellular automata can exhibit behavior that is far more complex than could be predicted from the simplicity of the individual cells. For a beautiful visual tour of cellular automata and a discussion of their actual and potential applications, see the book by Wolfram [97].