GPU Enables Obsession with Fractals
By Cleve Moler, MathWorks
I am a fractals addict. And I now have a GPU that is enabling my addiction. GPUs were originally intended to speed up graphics, but MATLAB^{®} uses them to speed up computation. Fractals^{1} are graphics that require extensive computation. Perfect candidates for a GPU.
The GPU has renewed my interest in the Mandelbrot set itself and facilitated work with three fractal variants: the “Burning Ship,” the “tower of powers,” and the global convergence behavior of Newton’s method.
My GPU is an NVIDIA^{®} Titan V, housed in a separate peripheral enclosure that is bigger than the laptop and that provides separate power and cooling. It is roughly 300 times faster than the CPU for computing these fractals, and completely changes how I interact with my mandelbrot
program. I introduced gpuArray
into the program and can now use a grid resolution comparable to my screen resolution and iterate until fine details in the image become visible.
Visualizing the Mandelbrot Set
The first published picture of the Mandelbrot set was produced on a line printer in 1978 by Robert W. Brooks and Peter Matelsk.
The first color image appeared on the cover of Scientific American in 1985. This was about the time that computer graphical displays were becoming widely available.
Since then, the Mandelbrot set has stimulated deep research in mathematics and has been the basis for numerous graphics projects, hardware demos, and web pages.
Inside the Mandelbrot Set
What happens in Mandelbrot stays in Mandelbrot. The Mandelbrot involves a simple iteration with complex numbers, starting at an initial point \(z_0\). The Mandelbrot set is the region in the complex plane consisting of the values \(z_0\) for which the trajectories defined by
\[z_{k+1} = z_{k}^2 + z_{0}, k=1, 2, ...\]
remain bounded. That's it. That's the entire definition. It's amazing that such a simple definition can produce such fascinating complexity.
The color image shown above is a bird’s eye view of the Mandelbrot set. It does not have the resolution to show the richly detailed structure of the fringe just outside the boundary of the set. In fact, the Mandelbrot set has tiny filaments reaching into the fringe region, even though the fringe is hard to see in this view.
Figure 1 is a sketch of the geometry of the Mandelbrot set. The largest component is a heartshaped cardioid, bounded by a curve with the parametric equation
\[z = e^{it}/2  e^{2it}/4, \pi \le t \le \pi \]
To the left of the cardioid is a circular disc of radius ¼. The cardioid and the disc touch at the point marked by the red dot. More about that red dot later.
Surrounding these two components are infinitely many nearly circular discs, and discs upon discs, with everdecreasing diameters. The exact locations and shapes of the smallest discs can only be determined computationally.
It has recently been proved that the Mandelbrot set is mathematically connected, but the connected region is sometimes so thin that we cannot resolve it on a graphics screen or even compute it in a reasonable amount of time.
Computing Mandelbrot
The fascinating patterns that we associate with Mandelbrot come from the fringe region just outside the set.
Here is the function that is the core of the computation. The input is a complex scalar starting point z0
and a real scalar parameter that I call depth
. The output is a count kz
. If the iteration is terminated because z
is outside the disc of radius of two, then z
was destined to become infinite and the count kz
can be used as an index into a colormap. On the other hand, if z
survives depth
iterations, then it is declared to be within the Mandelbrot set.
function kz = mandelbrot(z0,depth) z = z0; kz = 0; while (z*conj(z) <= 4) & (kz <= depth) kz = kz + 1; z = z*z + z0; end end
Let’s generate a grid of points in the complex plane covering a square of width w
, centered at the complex point zc
.
grid = 512; s = w*(1:1/grid:1); [u,v] = meshgrid(s+real(zc),s+imag(zc));
Now comes the only difference between a program that runs on the CPU and one that runs on the GPU. To generate the grid on a CPU,
z0 = u + i*v;
For the GPU,
z0 = gpuArray(u + i*v);
Then, for either processor, the statement
kz = arrayfun(@mandelbrot,z0);
applies the scalar function mandelbrot
to all the elements of the array z0
. On the CPU, this is a vectorized double for
loop over the grid. On the GPU, the statement is broken down into hundreds of individual tasks that run in parallel.
Now let’s look at two extraordinary images derived from the Mandelbrot fringe.
On the Fringe
The region known as the Valley of the Seahorses lies between the central cardioid and the disc to its left. There are actually two valleys, one above and one below the real axis. They meet at the red dot that we saw in Figure 1. With a little imagination, you can picture small marine creatures living in Figure 2. As we shall see later, the Valley also contains buried π.
Because the Mandelbrot set is selfsimilar, it contains an infinite number of miniature Mandelbrots, each with the same shape as the big one. The one shown in Figure 3 has a magnification factor of 10^{10}.
Buried π
A remarkable result was discovered in 1991 by Dave Boll, then a graduate student at Colorado State University. Boll was investigating the behavior of the Mandelbrot iteration in the Valley of the Seahorses. The valleys become narrower as they approach the axis, which they meet at (3/4,0), the red dot in Figure 1.
We can repeat Boll’s computation on a small grid centered just off the axis at the point 3/4 + yi for a tiny value of y and imaginary unit i.
y = 1.0e07 zc = 3/4 + y*i
Make the grid just large enough to touch the axis.
width = 2*y grid = 4
Choose depth
to be inversely proportional to y, which makes it huge.
depth = 4/y
With these parameters, run the code above. It produces
kz = 40000000 40000000 40000000 40000000 40000000 40000000 40000000 40000000 40000000 40000000 40000000 40000000 31415926 40000000 40000000 40000000 40000000 20943951 40000000 40000000 40000000 40000000 15707963 40000000 40000000
Look at the iteration count in the center of the grid. See a familiar value?
This isn’t a fluke. A 2001 paper by Aaron Klebanoff, “π in the Mandelbrot Set,” analyzes a similar computation in the cusp at the front of the cardioid.
Next, a curious variant of the Mandelbrot iteration.
The Burning Ship
The Burning Ship comes from a strange iteration:
\[z_{k+1} = F(z_k) + z_{0}, k=1, 2, ...\]
where
\[F(z) = (Re(z) + i Im(z))^2\]
I say this is strange because the function \(F(z)\) is not analytic. I am interested in this iteration because of the uncanny similarities in the following pictures.
The initial domain, shown in Figure 4, is a square of width 3.5 centered at 0.50.5i. I’ve inserted an arrow pointing to the region of interest in the wake of the ship.
Zoom in on the ship’s wake by a factor of 500 to the point 1.861.002i. Apply the bone
colormap to make it appear cold instead of burning. Figure 5 shows the resulting fractal next to a 1915 photograph of Antarctic explorer Ernest Shackleton’s ship Endurance frozen in the ice in the Weddell Sea.
Tower of Powers
Start with any complex number, \(z\), and repeatedly exponentiate it.
\[z, z^z, z^{z^z}, . . .\]
We can express this as an iteration. Start with \(y_0 = 1\) and let
\[y_{k+1} = z^{y_k}\]
If \(z\) is too big, this iteration will blow up to infinity. For some \(z\), it will converge to a finite limit. For example, if \(z = \sqrt 2\), the \(y_k\) will converge to 2. The most interesting case is when \(y_k\) approaches a cycle. For example, if \(z\) is near 2.5+4\(i\), the cycle has length 7.
2.4684 + 4.0754i 0.6216 + 0.3634i 0.2603  0.0184i 1.4868 + 0.3613i 3.4877 + 6.1054i 7.7632e06  2.6617e06i 1.0000 + 0.0000i 2.4684 + 4.0755i
This cycle length is the basis for the “tower of powers” fractal. Figure 6 shows the overall fractal.
Zooming in by a factor of 10^{5} and changing the color map produces the image shown in Figure 7.
Newton’s Method
When the starting point of Newton's method is not close to a zero of the function, the global behavior can appear to be unpredictable. Contour plots of iteration counts to convergence from a region of starting points in the complex plane generate thoughtprovoking fractal images.
The iteration is familiar. Pick a function \(f(z)\) with derivative \(f '(z)\). Start at \(z_0\) and let
\[z_{k+1} = z_{k} – f(z_{k})/f '(z_{k})\]
This will eventually converge to a zero of \(f\). Count the number of iterations it takes to get close.
There are many functions (and colormaps) to choose from. My favorite cubic polynomial is
\[f(z) = z^3 – 2z – 5\]
Figure 8 shows the complex plane divided into three regions where the iteration converges to one of the three zeroes of the cubic. Between these regions are areas of intense fractal action, shown in black in the figure.
Figure 9 shows the global behavior of Newton’s method seeking the zeroes of
\[f(z) = tan(sin(z))  sin(tan(z))\]
The function has infinitely many zeroes and an unbounded first derivative.
The function for Figure 10 is
\[f(z) = z\;sin(1/z) \]
The most prominent blue regions surround the zeroes at ±1/π.
Many of the images described here were new to me—I’d never seen them before. Interactive experiments with the GPU made them possible.
^{1} The term fractal was coined by Benoit Mandelbrot (1924–2010), a Polish/French/American mathematician who spent most of his career at the IBM Watson Research Center in Yorktown Heights, N.Y.
Published 2019
References

Aaron Klebanoff, “π in the Mandelbrot Set,” Fractals, World Scientific Publishing Company, vol. 9, 2001.

Frank Hurley, 1915, The long, long night [the Endurance in the Antarctic winter darkness, trapped in the Weddell Sea, Shackleton expedition, 27 August 1915], National Library of Australia.