Stories about Fractal Plotting · · · and HowTo



Revised in Summer of 2019


Contents Preparations: Canvases, Fractal Coloring, Dynamical Systems and Iterations
The Divergence Scheme
The Mandelbrot Set and Its Complexity
The General Convergence Scheme
Mandelbrot Fractals
Julia Fractals and Julia Sets
Newton Fractals
Fractal Enhancement
3D Rendering of Fractal Scenery
Speaking loosely without using technical terms such as the HausdorffBesicovitch dimension, a fractal is an object that is selfsimilar, i.e., a large part of it contains smaller parts that resemble the large part in some way; see Figures 0.10.5 below. Mathematician Benoit Mandelbrot coined the term “fractal” in 1975 and created a branch of mathematics called fractal geometry seven years later. As an "IBM Fellow," he had access to some of the best computers available for his research at the time.
Although fractal geometry is a very young field in 5,000year old mathematics, scientists have already discovered that it is filled with fascinating facts and applications as our world has fractals everywhere. For example, trees, mountains, blood vessels, mycelium strands, stock market graphs, weather patterns, seismic rhythms, ECG signals, brain waves and the Internet turned out to be all fractals.
The idea of fractal was not particularly new for the computer age, however, as Georg Cantor introduced the prototypical "Cantor Set" in 1883 almost 100 years before Mandelbrot published his book entitled "The Fractal Geometry of Nature." During the early 20th century, Pierre Fatou and Gaston Julia laid the foundations for fractals generated by "dynamical systems." It was in 1980 when Mandelbrot showed the famous fractal similar to Figure 0.1 called the "Mandelbrot set" generated by a simple dynamical system and a computer. Almost immediately after that, the novelty and complexity of the Mandelbrot set reinvigorated the interests in fractals and stimulated mathematicians to develop further theories in fractal geometry.
On the other hand, "chaos," often associated with fractals, was basically born as a brand new subject in 1974 from biologist Robert May's computer simulations of population dynamics through the dynamical system called the "logistic equation"  although some chaosrelated phenomena such as the sensitivity to a small change in the initial condition or the "butterfly effect" had been observed earlier by some mathematicians and physicists. Like "fractal," the word "chaos" was used as a mathematical term for the first time in 1975 when the American Mathematical Monthly published "Period Three Implies Chaos" by T.Y. Li and James Yorke. The paper received a great sensation especially because there appeared very little difference between chaotic and random outcomes even though the former resulted from deterministic processes. It is now known that fractals and chaos are closely related and together they provide applications in not only sciences but also art.
Painted on a sphere
Through Google, we find numerous websites that display stunningly beautiful computergenerated fractal art images. It indicates that a large population not only appreciates the digital art form but also participates in the eyeopening creative activity. Written below is a guide on how to program a computer and plot popular types of fractals generated by simple dynamical systems. It is not a text on computer programming /coding. Instead it tells the general principles needed for fractal plotting without going into too many specifics that might confine programmers' thought processes. It assumes the reader's basic programming experience and encourages him/her to be creative and engage in frequent computer experiments based on the essentials.
Particularly exciting is the moment the fractal image generated by our personal program emerges in our computer screen because of its rather unpredictable nature. The intricate fractal patterns often change dramatically when we alter some of the input values slightly  perhaps influenced by the butterfly effect and unpredictability of chaos. The reader who may be merely intrigued by the general idea behind fractal plotting is encouraged to try it. Many of the images will stir our imaginations in the part of mathematics that is in fact quite deep and still filled with unknowns. Best of all, though, it is plain fun.
§ 1. Preparations
The reader who wants to try it out needs to know (a) fundamental algebra and geometry of complex numbers, (b) beginning calculus, and (c) basic computer programming in such language as C++, Delphi and Java.
(a) includes the practice of writing a complex number z as a point (x,y) in the xyplane as well as the standard algebraic expression z = x + yi and ability to do basic arithmetic of complex numbers such as multiplication and division. The complex plane means the set of all complex numbers z = (x,y) which coincides with the Cartesian xyplane. For each complex number z = (x,y), the absolute value of z means z = √(x^{2} + y^{2}) and it represents the distance of z from the origin 0 of the complex plane. More generally, if z and w are complex numbers, z  w represents the distance between the complex numbers. (b) includes the basic ideas about the
derivative of a function such as a
critical point where the derivative vanishes. Particularly important in (c) is a twodimensional (2D) array of numbers.
We now introduce several preliminary ideas needed for fractal plotting.
Canvases: At the beginning of our fractal plotting process, we always lay out a rectangular canvas on which we compute a fractal using complex numbers. Therefore, a canvas is naturally a subset of the complex plane in the familiar Cartesian coordinate system, and we assume it is bounded by the numbers xmin, xmax, ymin and ymax, as shown in Figure 1.1 below. The fractal "goldfish" is just an example that can be plotted as output on the canvas at the end of the process.
Of course, our computers cannot handle infinitely many complex numbers in the canvas, so we specify the image resolution of the canvas using two integers imax and jmax and decompose the canvas into imax x jmax miniature rectangles called pixels or "picture elements." For example, if imax = 1200 and jmax = 800 then the canvas has the image resolution of 1,200 x 800 pixels and the output fractal will be plotted on the canvas comprising 960,000 pixels. Hence, our canvas is also a rectangle in the pixel coordinate system as shown in Figure 1.1, where each pixel of the canvas is represented by integral pixel coordinates (i, j).
To ensure that the two rectangles in Figure 1.1 are essentially the same, let
(1.1) Δx = (xmax  xmin) / imax; Δy = (ymax  ymin) / jmax,
and devise the following conversion formulas between the Cartesian and pixel coordinates:
(1.2) x = xmin + i Δx; y = ymax  j Δy;
(1.3) i = round[(x  xmin) / Δx]; j = round[(ymax  y) / Δy].
(1.2) allows us to choose exactly imax x jmax complex numbers in the first rectangle of Figure 1.1 in such a way that exactly one complex number (x, y) corresponds to each pixel (i, j) in the second rectangle. Using the onetoone correspondence, we shall henceforth assume that our canvas comprises these imax x jmax complex numbers and identify each complex number (x, y) of the canvas with the corresponding pixel (i, j) in the pixel coordinates. As a result, a pixel often becomes a complex number and vice versa in our upcoming arguments. Obviously, a canvas with a higher image resolution better "approximates" the "true" canvas containing infinitely many complex numbers.
Remarks: (1) The inputs to our computer program are a dynamical system (which we will discuss later in the preliminary section) and numbers imin, imax, xmin, xmax, ymin, ymax, and the output is a fractal image. Another and probably an easier way is to input the center and horizontal radius of the rectangle in Cartesian coordinates and program the computer to calculate xmin, xmax, ymin, ymax in such a way that the two rectangles in Figure 1.1 have the same aspect ratio.
(2) A good fractal image often comes from a microscopically "small" rectangle in Cartesian coordinates and a "large" canvas in pixel coordinates. That is because we often zoom in on a microscopic rectangular neighborhood in the complex plane using a high resolution canvas. In this text, a "large" canvas means a canvas with large imax and imin. For example, a canvas with 6,000 x 4,000 pixels is a large canvas.
Digital Fractals: If R is a canvas, a digital fractal on R simply means a 2D array (or matrix) of nonnegative integers over R. Since each element of R is a pixel (i, j) as we have just seen, we denote each component of the digital fractal by d(i, j). Our fractal plotting process comprises two phases, computing a "good" digital fractal and coloring it to produce a fractal image. The first phase will be described in § 2, § 4 and § 7 as a main body of the article.
Fractal Coloring: Although the coloring part occupies a relatively small part of our fractal plotting process, it is like painting a newly built house and it can improve or ruin the entire effort. It can also make two or more entirely different final images from basically the same digital fractal as we can see in the example below. The rightmost image shows that the second image, "Flock of Owls," is given by rotating the first image, "Racing Elephants," 90^{o} left and applying a different coloring. In our website, we use a rotation and a horizontal or vertical reflection of an image freely without warning.
Figure 1.2. Different Images from the Same Digital Fractal
Since a digital fractal is a two dimensional array of nonnegative integers over a canvas, our current goal is to convert nonnegative integers into colors that are likely to produce attractive fractal images.
When I began fractal plotting in the 1980s using Apple III and its Pascal compiler, there were only two colors, black and white, available and the coloring job was quite simple (although the fractal computation was excruciatingly slow). Practically every fractal emerged on the computer screen had a zebra pattern but I was still genuinely excited. By comparison, today's computers are so powerful they provide us with whopping 256^{3} = 16,777,216 colors that are packed in the color cube shown below and it requires certain mathematical tools to take full advantage of the technological treasure.
The 3D color cube has R, G, Baxes just like the familiar 3D space given by x, y, zaxes and each color is represented by an ordered triple (R, G, B) of "bytes" or integers between 0 and 255. For example, red = (255, 0, 0), green = (0, 255, 0), and
yellow = (255, 255, 0) = (255, 0, 0) + (0, 255, 0) = red + green.
Orange is the midpoint of the line segment joining red and yellow, namely (255, 128, 0). Also, black is the origin (0, 0, 0) of the color cube and white = (255, 255, 255). Thus a color is generally darker when all of its components are closer to 0. For example, a brownish color is given by darkening the components of orange like 0.7 x (255, 128, 0) = (179, 90, 0) and all shades of red are located on the Raxis between (0, 0, 0) and (255, 0, 0). We note that all of the color components are rounded to bytes between 0 and 255.
One of the most important objects in fractal coloring is a line segment joining two colors in the color cube: Suppose C_{1} = (R_{1}, G_{1}, B_{1}) and C_{2} = (R_{2}, G_{2}, B_{2}) are two colors in the color cube. Then the line segment from C_{1} to C_{2} is the set of all colors C = (R, G, B) satisfying the parametric equations
(1.4) R = (1  t) x R_{1} + t x R_{2}, G = (1  t) x G_{1} + t x G_{2},
B = (1  t) x B_{1} + t x B_{2},
where 0 < t < 1. Note that if parameter t = 0, C = C_{1}, and if t = 1, C = C_{2}, so we imagine that C moves along the line segment by changing its color gradually from C_{1} to C_{2} as "time" t progresses from 0 to 1.
To apply it on fractal coloring, suppose n_{1} and n_{2} are two nonnegative integers in the digital fractal such that n_{1} < n_{2} and we wish to color all numbers n between them. Let
(1.5) s = (n  n_{1})/(n_{2}  n_{1})
so s is between 0 and 1 for all pixels p on the canvas. Now set
(1.6) t = s
in (1.4) and color n accordingly. Then it is not difficult to check that as n changes from n_{1} to n_{2} the color placed on it changes from C_{1} to C_{2}. We also note that the graph of s vs n in (1.5) is a straight line segment slanting up from left to right as shown in Figure 1.4. We often modify the function s so its graph looks like the one shown below on the left:
Suppose C_{1} and C_{2} are yellow and brown, respectively. If we set t = s as in (1.6), then we get the color change that is shown in the figure on the left, and if we set t = s^{ 2} instead, then the "Vfunction" is converted to the parabolic "Ufunction" and we get a different color change as shown in the figure on the right. If we apply these two functions with a "correct frequency" on the digital fractal underlying Figure 1.4 above, we get the following "dragons," where the "background" is given by a straightforward color change from black to red.
If we use a cycle of blue → yellow → orange → red → blue → · · · instead of just two colors, then we get
corresponding to the functions
What I have shown above are just the basics of what I have developed on my own and people can extend them to their own coloring techniques. For example, the color cube has many other parametric curves in it and each one of them may provide an interesting coloring. One of the advantages in using computers is that we can find a lot of things such as a "correct frequency" just mentioned by a trialanderror process.
Once we get a satisfactory fractal, we may Photoshop to modify and/or improve it. I regularly use Photoshop's color and contrast editing functions. Actually, we can program our computer to do many of Photoshop's functions but cannot beat its quickness. I would say, however, that a poorly made fractal cannot be saved even by Photoshop, and we still have to learn how to get a good one in the first place.
Dynamical Systems and Iterations: When we solve a mathematical problem using a computer, we usually do it by exploiting what the machine does best, namely an iteration, which means repeating a certain process over and over, often for thousands or even millions of times, at a blinding speed. As an example of iterations, consider the equation called the Mandelbrot equation
(1.7)
z_{n+1} = z_{n}^{2} + p ,
which involves two indexed variables z_{n+1} and z_{n} and the third variable p called a parameter. All variables vary through complex numbers. The iteration index n is especially important for fractal plotting and it is there for us to iterate the equation to generate a sequence of numbers once the value of p and initial value z_{0} are given. For instance, let p = 2 and z_{0} = 0. Then setting the index n = 0, 1, 2, · · · in (1.7), our (properly programmed) computer iterates (1.7) and calculates the sequence of numbers
z_{1} = z_{0}^{2} + p = 0^{2}  2 = 2 ,
z_{2} = z_{1}^{2} + p = (2)^{2}  2 = 2 ,
z_{3} = z_{2}^{2} + p = 2^{2}  2 = 2 ,
and similarly, z_{n} = 2 for n = 4, 5, 6, · · · . If we hold the value of z_{0} at z_{0} = 0 and change the value of p from p = 2 to p = 1.9 in (1.7) then the computer again iterates (1.7) and calculates thousands of terms within a fraction of a second to give us the sequence of numbers
z_{0} = 0, z_{1} = 1.9, z_{2} = 1.71, z_{3} = 1.0241, · · · , z_{30} = 1.1626, z_{31} = 0.5483, · · · .
Similarly, if we leave the value of p fixed at p = 2 and change the value of z_{0} from z_{0} = 0 to z_{0} = 0.1, we get
z_{0} = 0.1, z_{1} = 1.99, z_{2} = 1.9601, z_{3} = 1.842, · · · , z_{30} = 0.7157, z_{31} = 1.4877, · · · .
Although we have shown only real sequences for simplicity, actual numbers involved in fractal plotting are complex numbers in the complex plane. Thus, the Mandelbrot equation (1.7) comprises infinitely many sequences of complex numbers, one sequence
z_{n} = z_{0}, z_{1}, z_{2}, z_{3}, · · ·
for each choice of values of p and z_{0}, that hop around all over the complex plane. It is therefore quite appropriate to call the Mandelbrot equation a "dynamic mathematical system" or dynamical system. As we shall see in § 5, there are infinitely many dynamical systems including (1.7) and the logistic equation (5.3). The next preliminary section shows how dynamical systems and digital fractals are connected.
Orbits, pCanvases and zCanvases: Consider the dynamical system (1.7) with a fixed value of z_{0}, say z_{0} = 0. As we have seen, each value of p in (1.7) generates a sequence z_{n} of complex numbers, which we call the orbit of p (with the fixed value of z_{0}). "Sequence" and "orbit" generally have a subtle difference in mathematics but we won't distinguish them here for simplicity.
Let R be a canvas, interpret it as a set of pvalues for the dynamical system, and call it a pcanvas. Then each complex number p in the pcanvas R generates a unique orbit z_{n} of p (with the fixed value of z_{0}). Let M be a large integer, say 10,000, and for each orbit z_{n} of p emanating from the pcanvas R, let m be the least iteration index less than M for which certain behavioral change (such as z_{m} > 100) occurs in the orbit. If no such index exists, let m = M. Because p is a pixel (i, j) as well as a complex number of the pcanvas (see Canvases) and m depends on p, we may write
m = d(p) = d(i, j).
We have just stated a basic principle for constructing a digital fractal comprising nonnegative integers d(i, j) for all p = (i, j) in the pcanvas R.
If we find a good digital fractal on a good canvas, it is likely that we get a good fractal image, and we will be talking about such digital fractals and canvases often in the upcoming sections. Figure 0.5 shows examples of fractals generated by the logistic equations (5.3) on pcanvases, and here are examples by the Mandelbrot equation (1.7):
Figure 1.6. Fractals on pCanvases Generated by the Mandelbrot Equation
By the symmetry argument, we can also talk about an orbit of z_{0} given by the dynamical system with a fixed value of p and a zcanvas comprising values of z_{0}. For each z_{0} = (i, j) belonging to the zcanvas, we may use the same argument and define a nonnegative integer m = d(z_{0}) = d(i, j). The 2D array comprising all d(z_{0}) = d(i, j) is a digital fractal on the zcanvas.
In Figure 1.7 below, we have examples of fractals, called "Julia fractals" in § 6, painted on zcanvases. The first two were generated by the Mandelbrot equation and the third by the logistic equation. We will discuss fractals on pcanvases in § 2  § 5 and those on zcanvases in § 6 and § 7.
Figure 1.7. Julia Fractals on zCanvases
§ 2. The Divergence Scheme (The Escape Time Algorithm)
We begin our discussion with the function f_{p}(z) = z^{2} + p involving parameter p so we can write the Mandelbrot Equation (1.7) as
(2.1) z_{n+1 }= f_{p}(z_{n}) = z_{n}^{2} + p.
The derivative of the function f_{p} is f_{p}'(z) = 2 z so the critical point of f_{p} is z = 0. Throughout § 2, the dynamical system means the Mandelbrot equation (2.1) with the fixed initial value z_{0} = 0, which is the critical point. Hence, all orbits of p are associated with the critical point z_{0} = 0 as well. As we shall see in § 6, such orbits play a crucial role in fractal geometry and is called critical orbits.
We say that a sequence z_{n} of complex numbers diverges to ∞ if the real sequence z_{n} diverges to ∞, i.e., if z_{n} gets further away from the origin of the complex plane without bound as n increases. Our current goal is to plot various fractals on pcanvases using the notion of divergence of orbits of p. Here is a classical example of plotting the most famous fractal called the Mandelbrot set, which is painted white in Figure 2.1.
Example 1: Let R be the square pcanvas with center 0 and radius 2.1, so it is bounded by xmin = 2.1, xmax = 2.1, ymin = 2.1, and ymax = 2.1, and suppose it has the image resolution of 400 x 400 = 160,000 pixels. Paint the whole canvas white initially.
Let M = 500 and for every pixel, er complex number (see Canvases), p in the pcanvas, iterate (2.1) at most M times. Thus, M is the maximum number of iterations that is intended to "approximate" ∞ and also rescue the computer from getting trapped in an infinite loop. Here's a simple instruction for our computer:
If
z_{1} > 2 then color the pixel p black ,
else if z_{2} > 2 then color the pixel p red ,
else if z_{3} > 2 then color the pixel p black ,
· · ·
else if z_{M} > 2 then color the pixel p red .
It can be shown that
(†)
z_{m} > 2 for some m if and only if the orbit z_{n} of p diverges
to ∞ ;
hence the above scheme assigns a color to each pixel p in the pcanvas R according to how quickly the orbit z_{n} of p escapes from the circle of radius 2 before taking a long journey toward ∞.
This basic scheme can be streamlined in a computer program by using:
just two variables zold and znew instead of the array z_{0} , z_{1} , z_{2} , ... , z_{n} , ... while keeping track of the important iteration index n.
z_{n}^{2 }> 4 in place of z_{n}
> 2 to avoid the hidden square root operation.
The Mandelbrot Set: By definition, the famous Mandelbrot set, or the M set for short, is the set of all complex parameters p whose critical orbits do not diverge to ∞, or equivalently, the set of all parameters whose critical orbits stay within the circle of radius 2 forever. The white "snowman" of Figure 2.1, which is left untouched by the above redblack algorithm and retains the initial white canvas color, depicts an approximation of the Mandelbrot set given by replacing "forever" by "up until M = 500 ". Theoretically, such approximation gets better as M gets bigger, but for our small canvas of 400 x 400 pixels, the use of any M > 500 will just waste the computing time without improving the output image.
Choosing pCanvases: (†) implies that the Mandelbrot set stays within the circle of radius 2, which explains why we chose the pcanvas containing the circle by a small margin in Example 1. Other pcanvases for (2.1) that produce good fractals are around the boundary of the Mandelbrot set. It is extremely unlikely to get something interesting if we randomly choose a pcanvas.
Computing Digital Fractals: Choose a pcanvas and a large integer M. For each p = (i, j) in the pcanvas, let
(2.2) d(i, j) = the least iteration index m < M such that the orbit of p satisfies z_{m} > 2
if such m exits, and if not, let d(i, j) = M. Then the 2D array comprising all d(i, j) is a digital fractal on the pcanvas. In our computer program based on (2.2), it is important to have the dynamical system stop iterating for each p as soon as the first (hence least) iteration index m is found so as to avoid wasting the computing time. If M and the pcanvas are as in Example 1, then we get Figure 2.1 by properly coloring the digital fractal.
We often refer to the iteration index m = d(i, j) as time needed for the orbit of p = (i, j) to escape from the circle of radius 2, hence the term escape time algorithm for the algorithm based on (2.2). The "escape time" has double meanings, however, and it sometimes means "time" needed for an algorithm to "escape" from the iteration loop so as to avoid wasting the computing time. The latter interpretation is much broader and applies to many computer programs. Since all of the escaping orbits diverge to ∞ by (†), we call the algorithm based on (2.2) the divergence scheme for fractal plotting to avoid the minor confusion. In § 4, we will use the convergence scheme which is also an escape time algorithm in the latter sense.
Example 2: Figure 2.2 is a cropped and resized image from a fractal on a pcanvas with 6400 x 3200 pixels and with center (0.2820607, 0.011014375) and horizontal radius 0.0000011. The center of the pcanvas is outside of the Mandelbrot set but very near the boundary. The large canvas and M = 100,000 are used to allow a high resolution printout. Part of the Mandelbrot set shown in the figure is painted black and the exterior of the Mandelbrot set various shades of green, sky blue, blue and red. Note that the figure shows several replicas of the snowman painted black that look like small isolated islands.
Figure 2.2. Fractal by the Mandelbrot Equation and the Divergence Scheme
Remarks: (1) Computers use most of the fractal plotting time for computing a digital fractal and it may take days and even weeks when the pcanvas and M are large. So, we almost always save the digital fractal on a hard disk before coloring it. Computers generally do the coloring job very quickly, although we humans may struggle for hours to figure out satisfying coloring.
(2) (†) implies that the use of any radius greater than 2, say 100, in place of 2 is valid in the divergence scheme, and it rarely degrades the output picture. Such an observation becomes important when we experiment with different dynamical systems without knowing what numbers to use in place of 2.
(3) As we have seen, the Mandelbrot set is based on the critical point of the function f_{p}, but it is not a requirement in a creative activity. In fact, the use of certain noncritical points provides us with interesting deforming effects on the Mandelbrot set.
.
§ 3. The Mandelbrot Set and Its Complexity
The Mandelbrot set is famous for a reason and it turned out to be one of the most complex figures ever plotted on a plane. The following celebrated theorem guarantees that no figures on the plane are more convoluted than the Mandelbrot set.
Shishikura's Theorem: The fractal dimension of the boundary of the Mandelbrot set is the same as the dimension of the plane, namely 2.
See Wikipedia for "fractal dimension." It also follows quickly from the theorem that the boundary of the Mandelbrot set has infinite length while it bounds a small area; see Figure 2.1.
If we alter the values of xmin, etc. and apply the divergence scheme with a greater number of colors on (2.1) to zoom in on a small area near the boundary of the "snowman" in Figure 2.1, we may catch a glimpse into how complex the set really is. For example, Figure 3.1 shows cropped/resized images from a fractal on a pcanvas with 4,000 x 4,000 pixels, which is given by zooming in on a microscopic rectangular neighborhood of the point
p = (0.25000316374967, 0.00000000895902)
with radius ≈ 0.0000000000001 = 10^{13} (instead of 2.1 used for Figure 2.1).
It shows a replica of the "snowman" at the center but the microscopic rectangle actually contains infinitely many such replicas. With experience our eyes become sharper in spotting where many of the replicas hide and it helps us realize that the Mandelbrot set is indeed selfsimilar and a fractal in that sense. By the definition of the Mandelbrot set, the "colorful" area (which is painted in colors other than black) depicts outside or the complement of the Mandelbrot set. Here is another amazing theorem:
The DouadyHubbard Theorem: The Mandelbrot set is connected.
So, all of the replicas of the snowman are connected by an invisible network that belongs to the Mandelbrot set and they dictate the formation of the intricate and complex pattern shown in Figure 3.1. What is also amazing is that the Mandelbrot set is "simply connected," i.e., the connecting network has no loops. It suggests that the complement of the Mandelbrot set comprising all those colorful pieces is also connected, and it has also been established as a theorem. There are still open questions involving such notions as "local connectedness" that defy the challenges by genius mathematicians.
Plotting the boundary of the snowman always demands a higher maximum number M of iterations because of ∞ involved in the definition of the Mandelbrot set, and for the picture on the left in Figure 3.1, we used whopping M = 1,500,000. If we use M = 500,000 (still a large number) instead, the outline becomes blurry as shown in the picture on the right. It is astounding that the difference between the two pictures comes from 500,000 to 1,000,000 additional iterations of the Mandelbrot equation for every black pixel.
Digression: I have two Pentium XP machines with 4 gigabytes of RAM and use one of them to do tedious and timeconsuming jobs. Many of the jobs involve a large canvas with as many as 6,400 x 6,400 pixels and often require days or even weeks of computing time. I bought the workhorse used for $89 by net shopping. Unlike in the 1980s when weaker computers were precious and expensive, people currently recycle or even throw away millions of powerful computers each year in the U.S.
We note that the body of the white snowman in Figure 2.1 is bounded by a "cardioid," whose cusp (or sharp point) is at the complex number (0.25, 0). It is very near the center (0.25000316374967, 0.00000000895902) of the zoomed area shown by Figure 3.1. To get an intricate picture, we almost always try to zoom in on a small neighborhood near the boundary of the big snowman. The closer the neighborhood is to the boundary, the more complex the fractal becomes. The reason lies in the aforementioned Shishikura's theorem and can also be explained in terms of chaos: the orbits of parameters p near the boundary have diverse properties and interact in a complex manner. A slight change in the value of p in the area may result in an utterly different escape speed of its orbit that dictates the pixel color.
Figure 3.2. A Pattern Near the Boundary of the Mandelbrot Set
Errors in Computation: Aside from programming bugs and other human errors, fractal plotting entails two unavoidable errors, each of which contributes to a loss of mathematical precision in the output. One is the truncation error resulting from "truncating" the infinite process after the finite number of steps given by M and the other is the roundoff error caused by our "imperfect" computer that has to round almost every real number involved. In my program, M is usually between 500 and 100,000 but it occasionally gets as large as 1,500,000 as shown in Figure 3.1. Using larger M is better in theory as it reduces the truncation error but is worse in practice as it makes the computing time longer and at the same time causes more roundoff errors to propagate. Balancing the good and the bad to find an optimum number M is a difficult problem in computing. Although it is not an issue if our goal is in art, it is something we should keep in our mind. Computers are great tools for mathematical research but they mislead us from time to time.
§ 4. The General Convergence Scheme
We are not done yet with the complexity of the Mandelbrot set and still stay with it. The Mandelbrot set has become so illustrious, everybody with at least some interest in fractals knows its shape by heart: The main body is in the form of a "cardioid" with a bunch of circular disks attached to it, and to each of the disks another bunch of disks are attached. The pattern repeats indefinitely as if the cardioid has children, grandchildren, great grandchildren, etc. etc., but beyond that nobody knows exactly what's happening.
Suppose we try to plot the Mandelbrot set on a sufficiently large pcanvas. By the definition of the Mandelbrot set, the orbits of all of the pixels belonging to the cardioid and the disks descending from it are bounded and the divergence scheme leaves all of the pixels in a single color like white (§ 2) and black (§ 3). So, the question is: Is there a way to make the Mandelbrot set more colorful?
It turned out that if p belongs to the interior of the cardioid, then the orbit of p is a convergent sequence in the complex plane, and how fast z_{n} converges to a point can be detected by
(4.1)
z_{n+1}  z_{n} < ε
if it is used like z_{n} > 2 in the divergence scheme of § 2. Here, ε is a fixed small number like ε = 10^{ 12}. Thus, the convergence scheme for computing the digital fractal for the area is given by
d(i, j) = the least iteration index m < M such that the orbit of p = (i, j) satisfies z_{m+1}  z_{m} < ε
instead of (2.2). The first inset in Figure 4.1 (where the Mandelbrot set is rotated to fit in the web page better) shows the interior of the cardioid painted by the convergence scheme with a color change from red to black.
Figure 4.1. The Mandelbrot Set by the General Convergence Scheme
The convergence scheme comes from the fact that if a sequence z_{n} converges then for any ε > 0, there is an integer N such that n > N implies (4.1). Unfortunately, the converse of the statement is not necessarily true, and we have to use the convergence scheme with care as it can be fooled by some of the nonconvergent orbits, notably orbits that diverge to ∞ very slowly. See Remark (2) at the end of the section.
It also turned out that if p belongs to the interior of the largest circular disk of the Mandelbrot set, then the orbit of p "eventually" becomes periodic of period 2 and starts to oscillate between two points in the complex plane. Like in the convergent case, the speed at which the orbit becomes periodic can be detected by
(4.2)
z_{n+k}  z_{n} < ε
with k = 2. The first inset of Figure 4.1 is a result of applying the general convergence scheme based on (4.2) with period index k = 1, 2, 3, · · ·. The computer program detected, for example, that the period of the interiors of the second largest disks (there are two of them) is 3 and painted them with a color change from orange to black just like the largest disk. The main image of Figure 4.1 is painted by the divergence scheme and the general convergence scheme; the latter is used for the Mandelbrot set in conjunction with the enhancing technique of § 8 to make the smaller disks more visible. Needless to say, the divergence scheme is used for the background (complement) of the Mandelbrot set.
Remarks: (1) The convergence scheme is a special case (with k = 1) of the general convergence scheme.
(2) As we mentioned earlier, if a fractal involves orbits, some of which diverge to ∞ while others converge, then the convergence scheme may be fooled by some of the divergent orbits and treat them as convergent orbits. We can avoid the situation by first applying the divergence scheme with large M to put aside "all" pixels whose orbits diverge to ∞ before applying the convergence scheme.
(3) The main image in Figure 4.1 is a result of a computer program that applies the divergence scheme first and then the general convergence scheme with k = 1, 2, 3, ... in that order. In fact, the program detects the "periodicities" k of many of the circular disks in the Mandelbrot set that have interesting numerical patterns. The numbers are neatly given by Wikipedia's Periodicity Diagram. It will become a handy tool when we plot "Julia sets" in § 6.
(4) The weakness of the general convergence scheme we discussed in (2) is not entirely bad as we can sometimes take advantage of it and plot interesting fractals.
(5) Finally, what happens if the divergence scheme and the general convergence scheme both fail and leave a hole in the fractal? The "black snowman" of Figure 3.1 is an example of such hole and there is no way of computing any periodicity on it if such exists. We can still paint the hole, however, by applying the "Euclidean metric," which we discuss in § 8. Figure 4.2, in which the fractals are painted on spheres for a change, shows fractals with a hole on the left and without on the right.
Figure 4.2. The Mandelbrot Set (Right) Colored by the Euclidean Distance of § 8
§ 5. Mandelbrot Fractals
Since it was published in 1980, the Mandelbrot set became so popular that a great many digital artists, mathematicians and computer programmers have explored around it and shown their fractal images on a variety of objects including webpages, posters, book covers, Tshirts, and coffee mugs. Although the complexity of the Mandelbrot set shown on a twodimensional canvas is boundless and the hidden beauty inexhaustible, it has become quite a challenge to unearth new patterns from the Mandelbrot equation (2.1) using available computers and software. Consequently, creative work calls for a modification of the formula, and there are infinitely many formulas available for it. We call a fractal given by a dynamical system of the form
(5.1)
z_{n+1 }= f_{p}(z_{n})
a Mandelbrot fractal if it is painted on a pcanvas, where f_{p} is any nonlinear function containing parameter p. The initial value for the orbits is normally (but not always) a critical point of f_{p}, so the orbits used in our computer program are normally (but not always) critical orbits of p emanating from the pcanvas.
The rather comical example shown above is a Mandelbrot fractal given by iterating
(5.2)
z_{n+1 }= f_{p}(z_{n}) = z_{n}^{3} + z_{n} + p
with z_{0} = i / √3 which is a critical point of f_{p}. The tip of the arrowhead is at the origin (0, 0) of the pcanvas and the tiny isolated figure on the right is another snowman (who propelled the arrow); see the inset for an enlarged image. Interestingly, the snowmen pop up in many fractals. Note: The picture is actually given by rotating the Mandelbrot fractal 90^{o} counterclockwise to fit in the computer screen.
The Logistic Equation: In 1838 Pierre Verhulst introduced a differential equation called the "logistic equation," which became a widely used mathematical model for population growth. If we replace the derivative in the equation by its approximating difference quotient to make it more suitable for computer applications, we get an equation of the form (5.1) given by
(5.3)
z_{n+1 }= f_{p}(z_{n}) = p(1  z_{n}) z_{n }.
The equation is again called the logistic equation (or the logistic "map") and is equally applicable in the population dynamics when the variables are restricted to real numbers. In 1974 while conducting computer simulations of certain population changes, biologist Robert May discovered "very complicated orbits" of (5.3), many of which appeared to change randomly and showed the butterfly effect. He called the phenomenon "chaotic" and soon after that the concept of "chaos" became a sensation in mathematics. In 1993 a "chaotician" appeared in Steven Spielburg's hit movie, "Jurassic Park," tacitly suggesting a possibility of chaos in the controlled dinosaur populations.
So, we naturally plot Mandelbrot fractals of the dynamical system (5.3) by expanding its variables to complex numbers; see Figure 0.5. The source of chaotic orbits Robert May discovered is the real interval [α, 4] with α ≈ 3.57, which is a part of the horizontal "antenna" at the right edge of the "global" image (first image) of Figure 0.5. Figure 5.2 above shows what we might see if we zoom in on a microscopic neighborhood of a point in the interval. In this example, a noncritical point z_{0} = 0.1 is used. The logistic equation (5.3) which originated from the Verhulst population model happens to be a close cousin of the Mandelbrot equation (0.4) and conceals fractals similar to the ones found near the Mandelbrot set as we can see in Gallery I.
Incidentally, people familiar with multivariable calculus may find a fun project in mapping an image like Figure 5.2 on surfaces like the ones shown below. We can see more examples of this sort in Gallery II.
§ 6. Julia Fractals and Julia Sets
Recall that a Mandelbrot fractal is generated by a dynamical system of the form (5.1) with a fixed initial value z_{0} and painted on a pcanvas. A fractal is called a Julia fractal if it is given by a dynamical system of the same form with a fixed parameter value p = p_{0} and painted on a zcanvas instead. Each complex value in the zcanvas works as an initial value z_{0} that generates a unique orbit of z_{0}, and just as before the orbits provide us with the divergence and general convergence schemes for plotting the Julia fractal. For example, consider the Mandelbrot equation
(6.1) z_{n+1 }= f_{p}(z_{n}) = z_{n}^{2} + p
with a fixed parameter value p. It is important to note that when z_{0} = 0, which is the critical point of f_{p}, the orbit of z_{0} coincides with the critical orbit of p used in § 2 to construct the Mandelbrot set.
Figure 6.1. Julia Fractals by the Mandelbrot Equation
Figure 6.1 shows two Julia fractals, along with their respective parameter values, painted on zcanvases centered at the critical point z_{0} = 0. The first parameter value comes from outside of the Mandelbrot set and the second from inside, or to be precise, inside of a circular disk of the Mandelbrot set of periodicity 11; see Wikipedia's Periodicity Diagram. The Julia fractal on the left is given by the divergence scheme alone and the one on the right by the divergence scheme followed by the general convergence scheme (4.2) with period index k = 11. I call the fractal on the right "Hydra of Lerma with Eleven Heads."
By the Julia set of (6.1) with a fixed value of p or simply the Julia set of p, we mean the boundary of the set consisting of all initial values z_{0} whose orbits diverge to ∞. In practical terms, the Julia set is the boundary of the subset of the Julia fractal left untouched by the divergence scheme. Hence, the Julia set is similar to the Mandelbrot set we saw in § 2  § 3 and dictates the formation of the Julia fractal. It is among the most important objects in fractal geometry as it displays a variety of fractal and topological structures and is the source of chaotic orbits that behave unpredictably.
Figure 6.1 shows Julia sets of nearly equal parameters with completely different topological structures. The boundary of the "Hydra" shown at right is the Julia set in the Julia fractal, and according to the upcoming theorem, it is connected as a set in one piece. On the other hand, the Julia set in the Julia fractal on the left happens to comprise "totally disconnected" powders and is called Cantor dust in topology. Georg Cantor, the pioneer of set theory, was among the first mathematicians who studied fractal structures in the 19th century.
The concept of Julia set appeared well before the computer era and it naturally extends to a more general function f_{p} like the idea of the Mandelbrot set. The following theorem established around 1920 explains why the Mandelbrot set is based on the critical orbits (and not on some other orbits) of p in the pcanvas.
The FatouJulia Theorem: Consider a dynamical system of the form
(6.2)
z_{n+1 }= f_{p}(z_{n}) = c_{n} z^{ n} + c_{n1} z^{ n1} + · · · + c_{2} z^{ 2} + c_{1} z + p ,
where c_{n}, c_{n1}, · · ·, c_{2}, c_{1} are constants. Then the Julia set of (6.2) with a fixed value of p is connected if and only if every critical orbit of p stays within a finite bound.
Note that if n > 2 in (6.2), p may have multiple critical orbits as f_{p} may have multiple critical points. In case of the Mandelbrot equation (6.1), each p has a single unique critical orbit; hence the FatouJulia theorem is beautifully stated as:
The Julia set of p is connected if and only if p belongs to the Mandelbrot set.
Furthermore, for the Mandelbrot equation we have:
If the Julia set of p is disconnected then it must be Cantor dust that is totally disconnected.
The four Julia fractals in Figures 6.2 and 6.3 below contain the Julia sets of (6.1) with parameters p taken from the interior of the Mandelbrot set. So, the Julia sets that bound the gold "lions" are all connected according to the FatouJulia Theorem, and consequently, the lions themselves are connected. (As mentioned earlier, we use rotations and horizontal/vertical reflections of fractal images freely.)
.
Figure 6.2. "Julia Lions" of Periods 14 and 42
Here's a popular way of plotting interesting Julia sets: Learn a simple pattern of the numbers in Wikipedia's Periodicity Diagram. Pick a parameter p in the interior of a disk of known periodicity and plot a Julia fractal of p on a zcanvas centered at the critical point 0. If the periodicity of the disk is k, paint inside the Julia set using the general convergence scheme (4.2) with period index k.
The left lion in Figure 6.2 is given by p = (0.296498, 0.020525) which is in a disk with k = 14 and the right lion by p = (0.296555, 0.020525) and k = 42. The two pvalues are very close to each other but they belong to different and adjacent disks. We can clearly see number 14 reflected in the shape of the left lion, but how about 42 = 14 x 3 in the right lion? People who find the lions in Figure 6.3 below share the same period 85 = 17 x 5 and spot that they are not mirror images of each other have the right idea and sharp eyes.
Figure 6.3. "Julia Lions" of What Period?
Additional examples of Julia sets we can find through the Mandelbrot set are shown in Gallery I. So, in addition to the astounding facts about the Mandelbrot set shown in § 3  § 4, we have just seen another: The Mandelbrot set conceals infinitely many varieties of Julia sets!
Figure 0.2 shown at the outset of this website is a Julia fractal generated by the dynamical system (5.2), for which the FatouJulia Theorem is applicable. Here, p = (0.185, 0.00007666), which is near the tip of the big arrowhead in Figure 5.1, and painted on a zcanvas containing the critical point z_{0} = i / √3 of the function. The gold "twin dragons" is given by the general convergence scheme (4.2) with period index k = 2 and the red background by the divergence scheme. If we alter the value of p and figure out the correct period, we get a variety of "twin dragons" including the pictures shown below. In general, the images are symmetric with respect to the center vertical line if p is a real number.
Figure 6.4. "Twin Dragons" of Dynamical System (5.2)
§ 7. Newton Fractals
A Julia fractal is called a Newton fractal if it is given by a dynamical system of the form
(7.1)
z_{n+1 }= z_{n}  g(z_{n})/g'(z_{n})
where p = 0 is invisible and g is an elementary function with its derivative g'. Although g is a function of a complex variable, the familiar rules of differentiation in high school calculus hold for g. In my program, g is almost always a polynomial which allows me to take advantage of the timesaving scheme called Horner's method (see the link at the end) to efficiently evaluate both g and g' that appear in the dynamical system. Horner's method is nothing but "synthetic division" taught in high school algebra, and it should be interesting for the reader to see how (differently) it is applied in computer programming.
The reader may have noted already that the dynamical system (7.1) is nothing but Newton's RootFinding Algorithm or Newton's Method for short, and therefore, each orbit converges to a root of the polynomial g quickly more often than it converges slowly or diverges. For this reason, we can plot most of the Newton fractals by the convergence scheme (4.1) alone with a relatively small maximum number of iterations M like 50  500. When we plot a Newton fractal, however, it is important that we take advantage of the roots of g to add more colors. The roots are what make the Newton fractal richer despite using only the convergence scheme with small M.
Coloring Newton Fractals: Suppose we already know the roots r_{1}, r_{2}, ..., r_{d} of the polynomial g, and choose mutually disjoint disks D_{1}, D_{2}, ..., D_{d} with centers r_{1}, r_{2}, ..., r_{d}, respectively. Initially, color the zcanvas say black.
Step 1. For each z_{0} in the zcanvas, use the convergence scheme (4.1) to iterate (7.1) and find the least iteration index m ≤ M such that the orbit z_{n} of z_{0} satisfies z_{m+1}  z_{m} < ε. Since the orbit is likely to converge, chances are that such m exists and z_{m} is in the disk about one of the roots. Go to Step 2 if m exists; do nothing otherwise.
Step 2. If z_{m} is inside D_{1} then paint the pixel z_{0} in a shade of say red and if it is in D_{2} then paint the pixel in a shade of say blue instead. Use the value of m to determine the amount of shading.
Figure 7.1. Newton Fractals of g(z) = z^{ 5}  1
"Crab Queue"
One of the simplest fractals for the beginners to draw is a Newton fractal of a polynomial of the form
g(z) = z^{ d}  1
as its roots, namely the dth roots of unity, are readily available.
The leftmost image of Figure 7.1 shows a Newton fractal of the polynomial with d = 5 given by the aforementioned coloring method. The zcanvas comprises 2,000 x 2,000 = 4,000,000 pixels centered at the origin with radius 1.1 and the maximum number of iterations is meager M = 100. Because of the small number of iterations for each orbit, my XP machine completed the fractal in a few minutes in spite of the fairly large canvas size. The second image that contains numerous "crab queues" is a slight variation of the first image.
The leftmost fractal consists of five regions and they are colored by various shades of sky blue, purple, red, amber, and blue. For example, the sky blue region comprises the initial values z_{0} whose orbits converge or "are attracted" to the root r = 1. The region is called the basin of attraction of Newton's method for the root r = 1. Thus, there are five basins of attraction in this example and their boundaries together form the basin boundary of Newton's method. The intricate boundary that separates the five regions shown in the leftmost fractal is actually a computer approximation of the basin boundary. The basin boundary is precisely the Julia set of the Newton fractal (which is a Julia fractal in the first place), and that is where Newton's rootfinding algorithm behaves in a "chaotic" fashion. Unlike the Julia sets we saw earlier, the Julia set of the Newton fractal is unbounded.
Remarks: (1) Some of the roots r_{1}, r_{2}, ..., r_{d} we discussed in the aforementioned coloring method may be repeated roots; hence the number of basins of attraction may be less than the number of roots. The same thing may happen, however, as a result of a programming error such as using disks D_{1}, D_{2}, ..., D_{d} that are not mutually disjoint. One way of avoiding the problem is to develop an alternative convergence scheme based on z_{n}  r_{k} < ε with k = 1, 2, ..., d instead of z_{n+1}  z_{n} < ε of the convergence scheme of § 4.
(2) If the initial canvas color is black as in the above example, a Newton fractal may contain "black holes" bounded by parts of the Julia set. It typically happens when some of the orbits become periodic of unknown periods. As illustrated in Figure 4.2,
one way of coloring such holes is to use the "Euclidean distance."
Figure 7.2. Cyclotomic Polynomials with Eight Roots
Another interesting example with known roots is a cyclotomic polynomial; see Wikipedia's
Cyclotomic Polynomial. The picture on the left in Figure 7.2 is a Newton fractal of the "30th cyclotomic polynomial"
g(z) = z^{ 8} + z^{ 7}  z^{ 5}  z^{ 4}  z^{ 3} + z + 1
with the unit disk highlighted. Since g happens to be a factor of z^{ 30}  1, its roots are among the 30th roots of unity that lie on the unit circle. In the picture, the thirty dots on the unit circle show where the roots of unity are located and eight of them colored yellow show the whereabouts of the roots of g. The picture on the right is a Newton fractal of the "20th cyclotomic polynomial"
g(z) = z^{ 8}  z^{ 6} + z^{ 4}  z^{ 2} + 1.
Cayley's Project: Figure 0.3 is another example and shows a Newton fractal of the cubic polynomial
g(z) = z^{ 3} + 1
painted on the unit disk. It has three basins of attraction as g has three distinct roots r_{1} = 1, r_{2} = (1  i√3)/2 and r_{3} = (1 + i√3)/2. The yellow basin, for example, comprises all initial values z_{0} in the zcanvas whose orbits z_{n} are attracted to r_{1} located at the leftmost point of the unit disk.
Arthur Cayley, one of the great mathematicians of the 19^{th} century who devised the idea of the basins of attraction of Newton's method, was frustrated by a problem he himself posed, which turned out to be equivalent to plotting the threecolor Newton fractal without a computer. The episode shows how lucky we are to live in the computer era as we have a huge advantage over the people born earlier in acquiring knowledge on complex structures. Even with the computer power, however, Newton's method still poses difficult basinrelated problems because of the unbounded and convoluted Julia sets. One such problem is to systematically find the roots of a polynomial by Newton's method if we don't know what they are.
Figure 7.3. Convoluted Newton Fractal on Plane, Sphere and Torus
Finding the Roots: As we have seen, a Newton fractal of a polynomial requires its roots. So, if we don't know them, what can we do? One way of finding them is to use Müller's algorithm; see e.g., Wikipedia's Müller's Method. Although Müller's method lacks the impressive simplicity and speed of Newton's rootfinding algorithm, it generally works well and automatically finds all roots of the polynomial. All Newton fractals shown in Gallery I use Müller's method even when the polynomials have known roots (so we don't have to list the roots in the computer programs).
For example, Figure 7.3 is a Newton fractal of a fifth degree polynomial whose roots are given by Müller's method. People who decide not to use Müller's method can still find the roots by Newton's method, although it may involve occasional "hitandmiss" routines. An important idea that applies to either method is to "deflate" the polynomial by synthetic division (Horner's method) each time a root is found so as to avoid catching the same root repeatedly. Incidentally, the last two pictures illustrate fun mapping projects for people familiar with multivariable calculus.
§ 8. Fractal Enhancement by Metrics, SemiMetrics and Directional Derivatives
Suppose a fractal we plotted has a "black hole" like the second image of Figure 8.1 below or it is somewhat flat and unappealing like image A in Figure 8.2. In this section, we'll develop several ways to improve such images, and it does not matter if they are on a pcanvas or a zcanvas.
Suppose p = (i, j) and q = (m, n) are pixels on a canvas R. Then by the Pythagorean Theorem, the distance between p and q is given by
(8.1) D(p, q) = p  q = √ [(i  m)^{ 2} + (j  n)^{ 2}].
If S is a nonempty subset of the canvas R, then for each pixel p on R, we define the distance from p to S by
(8.2) D(p, S) = minimum of D(p, q) where q is in S.
For each p = (i, j) on the canvas, let
(8.3) d(i, j) = round [D(p, S)].
This gives rise to a twodimensional array of nonnegative integers comprising d(i, j) on the canvas.
Example 1: Let S be the boundary of the Mandelbrot set in the second fractal shown below and compute a digital fractal d(i, j) by following (8.1)  (8.3). Then the image on the left is given by coloring the digital fractal in shades of green and sky blue. The numerous "fireflies" shown in the fractal are caused by the small (and nearly invisible) replicas of the "snowman" that are part of the computergenerated Mandelbrot set.
Figure 8.1. Application of Distance Function in Fractal Plotting
Partial Derivatives and Fractal Shading: Suppose R, S and D(p, q) are as in (8.1) and (8.2) and define d(i, j) by (8.3) without "round." Then assuming that i and j are continuous variables like real variables x and y, we can easily define the "partial derivatives" of d(i, j) as limits just like an ordinary derivative in high school calculus. What is important for our purpose is that, when i and j are integers, we can approximate the partial derivatives by the following simple formulas:
(8.9) P_{i}d(i, j) = d(i+1, j)  d(i, j) and P_{j}d(i, j) = d(i, j+1)  d(i, j).
For simplicity, we omit saying "approximation" and "derivative" and call P_{i}d and P_{j}d the partial of d with respect to (wrt) i and the partial of d wrt j, respectively. We can also interpret the partial derivatives just like the ordinary derivative. For example, P_{i}d(i, j) represents the rate of change of distance D of p from S at the pixel (i, j) per pixel length in "east" direction (or the direction of iaxis); see the pixel coordinate system in Figure 1.1. In particular, the partial P_{i}d(i, j) at the pixel (i, j) is positive (negative) if D increases (decreases) in "east" direction on a small neighborhood containing (i, j).
Example 2: Consider image A of Figure 8.2 below, which is a Julia fractal called "Hydra of Lerma with Eleven Heads." Our goal is to make the hydra look more threedimensional using the partial P_{j}d(i, j) pointing "south" direction (or the direction of jaxis).
Step 1: Use the digital fractal for image A to plot a blackandwhite image B, which has the welldefined boundary of the hydra.
Step 2: Let S be the boundary of the hydra and compute d(i, j) on the interior of S using (8.1), (8.2) and (8.3) without "round."
Step 3: Compute the partial P_{j}d(i, j) on the interior and boundary of S using (8.9).
Step 4: Normalize the values of P_{j}d(i, j) to the range between 0 and 1. Set P_{j}d(i, j) = 1 on the exterior of S.
Step 5: Shade image A by the partial derivative, i.e., multiply the color of image A at each pixel (i, j) on the canvas by the normalized value of P_{j}d(i, j) to get image D of Figure 8.2; see Fractal Coloring of § 1.
Remark: If we use a plain white canvas instead of image A in Step 5, then instead of image D, we get image C, namely the visualization of the partial derivative D_{j}k(i, j). I think it is itself an attractive fractal.
Figure 8.2. Fractal Shading by Partial Derivative
Directional Derivatives: If we examine image C of Figure 8.2 carefully, we can see the shading effect by the partial derivative in south direction. What if we want to use more general directions such as northwest to see different effects? Fortunately, we can easily modify the simple formulas in (8.9) to get a generalized partial derivative, called a directional derivative, representing the rate of change of D in a direction such as northwest. For example, the second Julia fractal of Figure 6.1 is shaded by the directional derivative of distance D from the boundary in southeast direction.
Metric and SemiMetric: There are numerous functions that behave like the distance (8.1) and that can be applied on fractal plotting with various effects. We can also use them in the next section when we discuss "3D rendering of fractal scenery."
A function D of two pixel variables is called a metric on a canvas R if it satisfies the following four conditions for all pixels p, q and r belonging to R:
(8.4) D(p, q) ≥ 0 ;
(8.5) D(p, q) = 0 if and only if p = q ;
(8.6) D(p, q) = D(q, p) ;
(8.7) D(p, q) ≤ D(p, r) + D(r, q) .
D is called a semimetric if it satisfies (8.4), (8.5), (8.6) but not necessarily the "triangle inequality" (8.7).
For example, the function (8.1) is a metric called the Euclidean metric as it coincides with Euclid's notion of distance. There are many other interesting metrics such as the taxicab metric, which is the distance conceived by taxicab drivers of Manhattan; see e.g., Wikipedia's Metric. Also, if c is a positive constant and D is a metric then cD is a metric. In the rightmost image of Figure 8.1, we colored the interior of the Mandelbrot set by using the digital fractal d(i, j) and its directional derivative in southeast direction based on the metric 10D where D is the Euclidean metric (8.1).
The Euclidean metric (8.1) without "square root," namely
(8.8) D(p, q) = (i  m)^{ 2} + (j  n)^{ 2},
is a semimetric but not a metric as it fails to satisfy the triangle inequality, and it can also be used for fractal enhancement and 3D rendering. The biggest advantage of (8.8) over (8.1) is that the computation of D(p, S) in (8.2) is so much faster and it can be used as a tool to develop useful digital fractals of the form
(8.9) d(i, j) = round [f(D(p, S))],
where f is a function of a single variable such as f(t) = 100√ t. Here's another "before/after" example of applying a semimetric and a directional derivative.
§ 9. 3D Rendering of Fractal Scenery
This section is for people who are familiar with the 3D xyzspace and have some experience with plotting the graph of a function like z = x^{2} + y^{2} by a computer. We recall that in such graph, z represents the altitude of the graph at (x, y). Since our discussions are based on "readymade" 2D fractals like in the previous section, we continue to use pixel coordinates (i, j) instead of (x, y). We will also use letter w for the third dimension of the 3D space and talk about the ijwspace instead of the xyzspace.
Given a 2D fractal plotted on a canvas R, our goal in the section is to plot a graph in the ijwspace which reflects some of the characteristics of the fractal, and we accomplish it by finding an appropriate altitude w of the graph at each pixel (i, j) on R.
The Basic Method: Start with a Mandelbrot, Julia or Newton fractal and suppose its digital fractal is given by d(i, j). Set
(9.1) w = d(i, j)
at each pixel (i, j) on the canvas R. This gives us a quick graph in the 3D space, but because of the chaosrelated nature of the fractal, the 3D graph tends to get literally "chaotic" with countless and abrupt ups and downs. The chaotic fluctuation of d(i, j)) generally results in intricate color changes in a 2D fractal but it is not welcomed in a 3D rendering. We can still find a few ad hoc ways to deal with them, however, and produce interesting fractal scenery.
Figure 9.1. "Mirage" in 3D Rendering
The desert scene on the left in Figure 9.1 is based on a Mandelbrot fractal given by the logistic equation and the divergence scheme. It "shows" several hills, but if we look at the picture on the right, we note that the hills are actually holes shown upside down. The undesirable chaotic part of the picture was conveniently thrown into the holes. Gallery I shows a few more examples given by the basic method.
The Basic Method with Interpolation: For people who know interpolation on a 2D regular grid of numbers such as "bicubic spline interpolation" and "bilinear interpolation," there are many possibilities in converting the 2D fractal into a 3D fractal. Applying such interpolation over the smaller number of d(i, j) chosen at the grid points (i, j) is an efficient way of taming the chaotic behavior of the function d.
Figure 9.2. Comparison of 2D and 3D Fractals
The image on the left in Figure 9.2 is a copy of "Logistic Elephants" and the one on the right is given by applying bicubic spline interpolation on its digital fractal at selected grid points. We can see, for example, that the "elephants' eyes" correspond to the steep hills in the 3D image. The 2D fractal has "infinitely many" such eyes but most of them are suppressed by the interpolation. The winter landscape of Figure 9.3 is similarly done on a Mandelbrot fractal.
Figure 9.3. "Snow Country" by Bucubic Spline Interpolation
The Popular Method: Perhaps the most commonly used 3D rendering is by means of metric and semimetric we discussed in § 8. For example, start with the Mandelbrot fractal of Figure 3.1, and let d(i, j) be the distance (Euclidean metric) of (i, j) from the boundary of the Mandelbrot set; see (8.1)  (8.3). Just as in (9.1), let
(9.2) w = d(i, j)
so, the altitude of the 3D graph at each pixel equals the distance of the pixel from the boundary of the Mandelbrot set. Then image A of Figure 9.4 is given by coloring the graph (9.2). The conical pits in the image are caused by the small replicas of "snowman" hidden in the Mandelbrot fractal we mentioned in § 3.
Image B is given by looking at image A from underneath and shows that the pits in image A are indeed conical. Image C is a combination of part of image A and part of image B with its viewing point at the far side of image B. Image D is based on the Mandelbrot set of Figure 8.1 and one of the semimetrics mentioned at the end of § 8. These semimetrics produce mountains with varying concavities.
Figure 9.4. 3D Rendering by Altitude = Function of Distance from the Boundary of the Mandelbrot Set
Figure 9.5 shown below is basically the same as image A of Figure 9.4 with minor changes in the metric and colors. All images of Figures 9.4 and 9.5, use directional derivatives to define sunny and shady sides. Probably the hardest part in plotting them is casting shadows of mountains on lower grounds. It does not use anything from higher mathematics but requires the persistence of the programmer.
Figure 9.5. "Mandelbrot Island"
Figure 9.6 shown below is similar to image B of Figure 9.4 and based on the Mandelbrot set of Figure 8.1. Its colors are directly transferred from the second image of Figure 8.1. We again note that the numerous mountain peaks of the 3D image show exactly where the small replicas of the snowman hide in the original 2D fractal.
Figure 9.6. "Mandelbrot Volcano"
Finally, how about using metrics that are unrelated to the Euclidean distance? The picture shown below at left uses the "maximum metric" and the one at right the "taxicab metric" mentioned in § 8. The "Crater Lake" on top of the "Taxicab Mountain" is given by the Basic Method applied on the Mandelbrot set of Figure 0.1.
