Stories about Fractal Plotting · · · and HowTo




Contents Preparations: Dynamical Systems and Iterations, Canvases, and Fractal Coloring
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 a geometric figure that is selfsimilar, i.e., a large part of it contains a great many smaller parts that resemble the large part in some way; see Figures 0.10.5 below. Nature is filled with fractals as seen in such objects as mountains, shorelines, trees (branches, barks and roots), ferns, fluid flow patterns, cloud formations, blood vessels and mycelium strands. Benoit Mandelbrot coined the term "fractal" in 1975 and published a book entitled "The Fractal Geometry of Nature" in 1982. As an "IBM Fellow," he had access to some of the best computers available at the time for his research.
The idea of fractal was not particularly new for the computer age, however, as Georg Cantor introduced the prototypical "Cantor Set" in 1883  almost exactly a hundred years prior to Mandelbrot's book. 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 surprisingly many applications in sciences  and 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 the 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 array of numbers.
We now introduce several preliminary ideas needed for fractal plotting.
Dynamical Systems and Iterations: When we solve a mathematical problem using a computer, we often do it by exploiting what the machine does best, namely iterations, which means repeating a certain process over and over for thousands or even millions of times, at blinding speed. For example, consider the equation called the Mandelbrot Equation
(1.1)
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.1), we get 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 ,
which shows z_{n} = 2 for all n > 1. If we fix the value of z_{0} at z_{0} = 0 and change the value of p from p = 2 to p = 1.9 then our computer 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.
Remark: Let R be a rectangular area in the complex plane. Then the above computation shows:
(a) We can iterate the equation using a fixed value of z_{0} and "every" complex number p varying through R ;
(b) We can iterate the equation using a fixed value of p and "every" complex number z_{0} varying through R .
The observations (a) and (b) will become important later. We also observe that the Mandelbrot equation (1.1) comprises infinitely many sequences, one sequence for each choice of values of p and z_{0}, that hop around all over the complex plane. So it is quite appropriate to call the Mandelbrot equation a "dynamic mathematical system" or dynamical system.
As we'll see in § 5, there are infinitely many dynamical systems including the logistic equation (5.3), and at the outset of our fractal plotting program, we always furnish one of them as it determines our output image. For example, the logistc equation produces, among many others, the fractals shown in Figure 0.5, and the endless list of fractals by the Mandelbrot equation (1.1) includes Figure 0.1 and the following images:
Figure 1.1. Fractals Generated by the Mandelbrot Equation
We note that each of the fractals is painted on a particular "canvas" by particular coloring. In the next two topics, we discuss "canvases" and "fractal coloring."
Canvases and the Conversion Formulas: At the beginning of our fractal plotting process, we not only furnish a dynamical system but also lay out two rectangular canvases like the ones shown in Figure 1.2 below. One is in the complex xyplane and the other in the pixel coordinates (i, j) in computer graphics. That's because the fractal plotting we consider consists of two phases, computing a precolored fractal (a fractal in digital form) and coloring the fractal. The computing part uses the dynamical system of complex numbers and it is naturally done on a canvas in the complex plane, while the coloring part is done by working with pixels (short for "picture elements") in pixel coordinates. Note that the two canvases in Figure 1.2 show identical fractal images, but imagine that the fractal on the Cartesian canvas is actually a precolored fractal comprising a bunch of nonnegative integers instead of colors.
Roughly speaking, we compute a precolored fractal by interpreting the Cartesian canvas as a set of parameter values p and iterating the dynamical system for each p in the canvas; see Remark (a). Then as we'll see in § 2 and others that follow, the precolored fractal is given by the iteration indices n for which z_{n} shows certain behavioral change. There is an optional process, in which we interpret the Cartesian canvas as a set of initial values z_{0}; see Remark (b). We call the Cartesian canvas as a set of parameters p a pcanvas and the canvas as a set of initial values z_{0} a zcanvas.
For simplicity, let's assume that the Cartesian canvas means say a zcanvas for now. The input numbers xmin, xmax, ymin and ymax shown in Figure 1.2 define the boundary of the Cartesian canvas and imax and jmax define the size of the pixel canvas. The greater the size, the better the resolution of the printout of the fractal image computed on the Cartesian canvas.
The following conversion formulas define the important relation between the two canvases and allow our computer to move from the aforementioned first phase to the second smoothly.
(1.2) Δx = (xmax  xmin) / imax; Δy = (ymax  ymin) / jmax;
(1.3) x = xmin + i Δx; y = ymax  j Δy;
(1.4) i = round[(x  xmin) / Δx]; j = round[(ymax  y) / Δy].
Because the pixel coordinates are all integers, the pixel canvas consists of only finitely many pixels (i, j), and as a necessary step in computer programming, formula (1.3) conveniently chooses finitely many complex numbers (x, y) in the Cartesian canvas that are in a onetoone correspondence with the pixels in the pixel canvas. These are the only complex numbers we consider in the Cartesian canvas, and hence, we assume that the two canvases comprise exactly the same number of elements.
When it comes to our discussion, we will henceforth
(◊) identify the pixel canvas with the Cartesian canvas and each pixel with the corresponding complex number,
for the sake of simplicity. Thus, we'll often use the words "pixel" and "complex number" interchangeably when we talk about the canvas.
Remark on Aspect Ratios: The conversion formulas do not impose any restrictions on the aspect ratios of the two canvases to provide flexibility, but we may want the ratios to be equal to get mathematically accurate output. If that's the case, we can input imax and jmax and the center and horizontal radius of the Cartesian canvas and instruct our computer to find appropriate xmin, xmax, ymin, ymax.
Fractal Coloring: As we mentioned earlier, our fractal plotting program has two phases, computing a precolored fractal and coloring it. The first step is described in § 2, § 4 and § 7, so let's assume we are done with it for now and discuss the second step, namely, coloring a precolored fractal. As we shall see, a precolored fractal is a two dimensional array (or matrix) of nonnegative integers over our canvas, so our current object is to convert nonnegative integers into colors that produce attractive colored fractals.
The coloring part 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 products from the same precolored fractal as we can see in the above example. 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 choice. In our website, we use a rotation and a horizontal or vertical reflection of an image without warning.
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) satifying the parametric equations
(1.5) 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 numbers in the precolored fractal such that n_{1} < n_{2} and we wish to color all numbers n between them. Let
(1.6) 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.7) t = s
in (1.5) 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.6) 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.7), 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 on the fractal in 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
§ 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.1) as
(2.1) z_{n+1 }= f_{p}(z_{n}) = z_{n}^{2} + p.
Our goal is to plot various fractals using (2.1). As we have seen, once a value of z_{0} is given in (2.1), each value of p determines a unique sequence z_{n} of complex numbers, which we call the orbit of p with initial value z_{0}. "Sequence" and "orbit" generally have a subtle difference in mathematics but we won't distinguish them here for simplicity. We say that an orbit z_{n} of p diverges to ∞ if the absolute value z_{n} diverges to ∞, i.e., if z_{n} gets further away from the origin of the complex plane without bound as n increases. In § 2, we examine the real sequence z_{n} to find a behavioral change of the orbit that contributes to the fractal computation.
The derivative of the function f_{p} is f_{p}'(z) = 2 z so its critical point is z = 0. As we shall see in § 6, the orbit of p with initial value z_{0} = 0, which is the critical point, plays a crucial role in fractal geometry, so let's call it the critical orbit of p. Set z_{0} = 0 so all orbits are critical orbits throughout § 2.
Example 1: Let R be the square pcanvas with center 0 and radius 2.5, so it is bounded by xmin = 2.5, xmax = 2.5, ymin = 2.5, and ymax = 2.5, and suppose it comprises 400 x 400 = 160,000 pixels, er complex numbers; see (◊) in the previous section. Paint the whole canvas white initially.
Let M = 500 and for every pixel 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 the computer to plot Figure 2.1 shown below:
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, the use of any M > 500 will just waste the computing time without improving the output image.
PreColored Fractals: For each pixel p = (x, y) in the pcanvas, let
(2.2) k(x, y) = least iteration index m ≤ M such that the orbit of p satisfies z_{m} > 2.
Then the twodimensional array comprising all k(x, y) is the precolored fractal for the Mandelbrot set of Figure 2.1. To get the colored fractal, all we need is to have the computer color p = (x, y) red or black depending on the parity of k(x, y). It is important to note that the algorithm based on (2.2) stops the iteration of (2.1) as soon as the first (hence least) iteration index m is found so as to avoid wasting the computing time.
We often refer to m = k(x, y) as time needed for the orbit of a parameter p = (x, y) 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 the 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 6400 x 3200 pcanvas with center (0.2820607, 0.011014375) and horizontal radius 0.0000011. The large canvas and M = 100,000 are used to allow a high resolution printout. The figure shows part of the Mandelbrot set given by the black canvas color left untouched by the divergence scheme. Here, various shades of green, sky blue, blue and red are used instead of red and black in Example 1.
Figure 2.2. Fractal by the Mandelbrot Equation and the Divergence Scheme
Remarks: (1) Computers use most of the fractal plotting time for computing the precolored fractal and it may take days and even weeks when the pcanvas and M are large. So, we almost always save the precolored fractal on a hard disk before trying out various colors on it. Computers generally do the coloring job very quickly, although we humans may struggle for hours to get a satisfying color combination.
(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) The playground of the orbit of p for each p in the pcanvas is the entire complex plane and not at all confined to the pcanvas.
(4) 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.
(For "fractal dimension," see Wikipedia.) 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 large 4,000 x 4,000 canvas, 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.5 used in Figure 2.1).
It shows a replica of the "snowman" at the center but the rectangle of the atomic size 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 therem. Most of the mathematicians in fractal geometry conjecture that the Mandelbrot set is also "locally connected," but it has so far defied the efforts by genius researchers. If it is true, then every point in the Mandelbrot set has a neighborhood in which the local portion of the Mandelbrot set is connected just like the global Mandelbrot set of Figure 2.1.
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, grandchildern, 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 the prescribed canvas 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^{ 8}. Thus, the convergence scheme for computing the precolored fractal for the area is given by
k(x, y) = least iteration index m ≤ M such that the orbit of p = (x, y) 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. For example, the "eyes" of the "cuttlefish" in Figure 6.5 are painted by the convergence scheme that are fooled by the orbits that diverge to ∞ slowly.
(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 a nonlinear function containing parameter p such as f_{p}(z) = z^{3} + z^{2} + p , f_{p}(z) = z^{6} + p , and f_{p}(z) = p(z^{5}  z) + 1. The initial value for the orbits is normally (but not always) a critical point of f_{p}.
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 rondomly 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 sequence z_{n} called the 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 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 detect 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 like the one on the left if p is a real number.
Figure 6.4. "Twin Dragons" of Dynamical System (5.2)
Relation Between zCanvas and pCanvas: We have seen that given a dynamical system like z_{n+1} = z_{n}^{2} + p , there are two ways of plotting its fractal: Choose a value of p, say p_{0}, and plot the fractal on a zcanvas or choose a value of z_{0} and plot it on a pcanvas. The former produces a Julia fractal and the latter a Mandelbrot fractal. We also recall that a pcanvas is a set of pvalues or the points (p.r, p.i) where p.r and p.i denote the real and imaginary parts of p, and likewise for a zcanvas.
Since a pcanvas is always associated with a fixed z_{0} = (z_{0}.r, z_{0},i), it is more convenient to write the pcanvas as a set of ordered quadruples (z_{0}.r, z_{0}.i, p.r, p.i), and similarly, zcanvas as a set of ordered quadruples (z.r, z.i, p_{0}.r, p_{0}.i). Because z_{0} and p_{0} are fixed values, it follows that the p and zcanvases are planes in the four dimensional space comprising all points of the form (z.r, z.i, p.r, p.i). Furthermore, people familiar with multivariable calculus find that the two canvases are orthogonal to each other in the space and have a point in common, namely, (z_{0}.r, z_{0}.i, p_{0}.r, p_{0}.i).
Example: Consider the Mandelbrot equation z_{n+1} = z_{n}^{2} + p and let z_{0} = (z_{0}.r, z_{0},i) = (0, 0) and p_{0} = (p_{0}.r, p_{0},i) = (0.25000316374967, 0.00000000895902). Figure 6.5 shows Mandelbrot and Julia fractals plotted on microscopic pcanvas (left) and zcanvas, respectively, with the same center (z_{0}.r, z_{0}.i, p_{0}.r, p_{0}.i).
Figure 6.5. "Partying Cuttlefish" around the Mandelbrot Set
§ 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 a fun project 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 formulate several ways to improve such an image. Since a canvas containing a fractal is in the pixel coordinates (i, j), it is now more convenient and natural that we write every complex number, er pixel, as (i, j) rather than (x, y) in the Cartesian coordinates; see Canvas and Conversion Formulas of § 1.
Suppose p = (i, j) and q = (m, n) are pixels on a canvas R. Then by the remark shown at the outset of § 1, 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.
In particular, if S is itself a fractal such as a Julia set in a Julia fractal on a zcanvas, let
(8.3) k(i, j) = round [d(p, S)]
for each p = (i, j) on the canvas. Then the twodimensional array comprising k(i, j) is a "new" precolored fractal given by the "old" fractal S and the distance function (8.1).
Example 1: Let S be the boundary of the Mandelbrot set in the second fractal shown below and compute the precolored fractal k(i, j) by following (8.1)  (8.3). Then the image on the left is given by coloring the precolored 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 k(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 k(i, j) as limits just like an ordinary derivative in high school calculus; e.g., see Math Insight. What is important for our purpose is that, when i and j are integers, we can approximate them by the following simple formulas:
(8.9) D_{i}k(i, j) = k(i+1, j)  k(i, j) and D_{j}k(i, j) = k(i, j+1)  k(i, j).
For simplicity, we omit "approximation" and sometimes even "derivative" and call D_{i}k(i, j) and D_{j}k(i, j) the partial of k with respect to (wrt) i and the partial of k wrt j, respectively. The partial derivatives can also be interpreted just like the ordinary derivative. For example, D_{i}k(i, j) represents the rate of change of distance k from S per pixel length in "east" direction (or the direction of iaxis); see the pixel coordinate system in Figure 1.2. Thus, the partial D_{i}k(i, j) is positive (negative) if k 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 D_{j}k(i, j) pointing "south" direction (or the direction of jaxis).
Step 1: Use the precolored fractal for image A to plot a blackandwhite image B, which has a welldefined boundary of the hydra.
Step 2: Let S be the boundary of the hydra and compute k(i, j) on the interior of S using (8.1), (8.2) and (8.3) without "round."
Step 3: Compute the partial D_{j}k(i, j) on the interior of S using (8.9).
Step 4: Use a formula similar to (1.6) to "normalize" the values of D_{j}k(i, j) to the range between 0 and 1. Set D_{j}k(i, j) = 1 on the complement of the interior of S.
Step 5: Multiply the color of image A at each pixel (i, j) on the canvas by D_{j}k(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 as the rate of change of distance from the boundary 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 k in a direction such as northwest. For example, the second Julia fractal of Figure 6.1 is shaded by the directional derivative of distance 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 in the way shown in Examples 1 and 2 with varying 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(p, q) is a metric then cd(p,q) is a metric. In the rightmost image of Figure 8.1, we colored the interior of the Mandelbrot set by using k(i, j) and its directional derivative in southeast direction based on 10d(p, q) where d(p, q) is the Euclidean metric.
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, but it can also be used for fractal enhancement and 3D rendering. We also note that computation of (8.2) involving "square root" and "minimum" may easily take days, but it can be considerably shortened by the use of (8.8). We simply drop "square root" from (8.1) and (8.2) and replace (8.3) by
(8.3)' k(i, j) = round [d(p, S)^{ 1/2}].
In fact, if we replace d(p, S)^{ 1/2} by d(p, S)^{ s} with positive constants s or more general functions f(d(p, S)) in (8.3)', we create a host of precolored fractals k(i, j) based on distinct semimetrics. 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 precolored fractal is given by k(i, j). Set
(9.1) w = k(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 k(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. I sometimes wonder if some of the computergenerated 3D movie scenery are based on illusions like our mirage.
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 k(i, j) chosen at the grid points (i, j) is an efficient way of taming the chaotic behavior of the function k.
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 precolored 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 k(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 = k(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 mountain peaks of the 3D image show exactly where the small replicas of 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 "Taxicab's Mountain" is given by the Basic Method applied on the Mandelbrot set of Figure 0.1.
