Stories about Fractal Plotting
Contents Preparations: Canvases, Dynamical Systems and Iterations
The Divergence Scheme
The Mandelbrot Set
The Convergence Schemes
Mandelbrot Fractals
Julia Fractals and Julia Sets
Newton Fractals
Supplement
Fractal Coloring
Galleries
Gallery 2D
Gallery 3D
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.6 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.
Our world has fractals everywhere exemplified by trees, mountains, blood vessels, mycelium strands, stock market graphs, weather patterns, seismic rhythms, ECG signals and brain waves. In its article entitled "How Mandelbrot's fractals changed the world," the BBC states that fractal geometry has practical applications in diverse areas including diagnosing some diseases, computer file compression systems and the architecture of the networks that make up the Internet.
Figure 0.1. The Mandelbrot Set
 Figure 0.2. Julia Set
 Figure 0.3. Newton Fractal

The idea of fractal was not particularly new in mathematics for Mandelbrot's time or the computer age, 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 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.
 Figure 0.4. Classical Fractals by Geometric Constructions


Koch Snowflake, 1904
 Sierpinski Rectangle, 1916
 Pythagorean Tree, 1942

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. 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 soon became known that fractals and chaos are closely related and together they provide applications in sciences as well as in art.
Figure 0.5. Fractal from the Area of Chaos
NOVA broadcast on PBS in 2008: Largely because of its haunting beauty, the Mandelbrot set has become the most famous object in modern mathematics. It is also the breeding ground for the world's most famous fractals. Since 1980, the set has provided an inspiration for artists, a source of wonder for schoolchildren, and a fertile testing ground for the science of linear dynamics.
Figure 0.6. From the Boundary of the Mandelbrot Set
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 or coding. Instead it tells the general principles needed for fractal plotting without going into too many specifics. It assumes the readers' basic programming experience and encourages them 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―influenced by the "butterfly effect" of chaos, perhaps. 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 readers who want to try it out need to know (a) fundamental algebra and geometry of complex numbers, (b) beginning calculus, and (c) basic computer programming in such languages as C++, Delphi (Pascal) 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, division and exponentiation. 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 O 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 and a critical point of a function where the derivative vanishes. Particularly important in (c) is a twodimensional (2D) array of real numbers, which will be used every time a fractal is generated.
We now introduce several preliminary ideas needed for fractal plotting.
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.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 complex 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), our properly programmed computer iterates (1.1) 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.1) then the computer again iterates (1.1) 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, · · · .
Note that the behavior of the sequence may change drastically if we alter the value of p or z_{0} slightly. We exploit such changes to draw an intricate fractal with a variety of colors.
We have shown only real sequences for simplicity, but actual numbers involved in fractal plotting are complex numbers in the complex plane. Thus, the Mandelbrot equation (1.1) comprises infinitely many sequences of complex numbers, one sequence z_{n} for each choice of values of p and z_{0}. Because most of the infinitely many sequences dance around in the complex plane as the iteration index n increases, it is 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.1) and the logistic equation (5.3), each of which generates infinitely many fractals.
Figure 1.1. Fractal Generated by the Logistic Equation
Canvases: We begin with a simple example.
Let R be the rectangle in the complex plane defined by 2 ≤ x ≤ 2 and 1.28 ≤ y ≤ 1.28 and suppose we wish to plot the graph of the inequality x^{2} + y^{2} ≤ 1 on R using a computer. We first decompose R into, say, 50 × 32 miniature rectangles of equal size called picture elements or pixels and then represent the pixels by pixel coordinates (i, j) in such a way that the upper left and lower right pixels are (0, 0) and (49, 31), respectively. Thus, the i and jaxes of the pixel coordinate system are the rays emanating from the upper left corner of R and pointing east and south, respectively; see the diagram in Figure 1.2 on the left.
Let imax = 50, jmax = 32, xmin = 2, xmax = 2, ymin = 1.28 and ymax = 1.28. Then for each i = 0, 1, 2, · · ·, imax1 and j = 0, 1, 2, · · ·, jmax1, the pixel (i, j), which is a rectangle, contains infinitely many complex numbers (x, y). For our computational purpose, we choose exactly one representative complex number (x, y) in the pixel (i, j) by setting
(1.2) Δx = (xmax  xmin) / imax; Δy = (ymax  ymin) / jmax,
(1.3) x = xmin + i Δx; y = ymax  j Δy.
Consequently, we may view R as the rectangle comprising imax x jmax = 50 x 32 pixels, each of which has a unique representative complex number. The rectangle R with the pixel structure is called a canvas for plotting the output image with the image resolution of 50 × 32 pixels.
Plotting the graph of the inequality x^{2} + y^{2} ≤ 1 on the canvas is now easy. We examine each pixel (i, j), er the representative complex number (x, y) on the canvas R. If it satisfies the inequality, color the pixel red, and otherwise, color it white. Since the coloring process uses only finitely many pixels of the canvas R, the output image that resembles the Japanese flag is an approximation of the true graph. The greater the number of pixels, the higher the image resolution and the more accurate the output image.
Figure 1.2 shows two approximations of a fractal called "Goldfish in Love." The one on the left is painted on a canvas with 50 × 32 pixels and the other on a canvas with 500 × 320 pixels.
Figure 1.2. Fractal with Different Image Resolutions
In general, we can define a canvas using any positive integers imax and jmax and any real numbers xmin, xmax, ymin and ymax satisfying xmin < xmax, ymin < ymax and (ymax  ymin)/(xmax  xmin) = jmax/imax. The last equality assures that the rectangle has the correct aspect ratio (so the red circle in the Japanese flag output would not look oval).
pCanvases, zCanvases, Orbits, and a Preview of § 2  § 7: As we have seen with the Mandelbrot equation (1.1), a dynamical system comprises infinitely many sequences of complex numbers, one for each choice of values of z_{0} and p. For example, if we fix the value of z_{0} at, say, z_{0} = 0, then each value of the parameter p gives rise to 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.
We have also seen that a canvas is a rectangle in the pixel coordinate system comprising imax × jmax pixels for some positive integers imax and jmax and that each pixel on the canvas is identified with a unique representative complex number belonging to the pixel. Interpreting the complex numbers representing the pixels as values of p, we call the canvas a pcanvas for the dynamical system (with the fixed value of z_{0}).
Each pixel, er parameter p, on the pcanvas corresponds to a unique orbit of p that usually dances around in the complex plane, and as we shall see in § 2  § 5, we plot a fractal called a "Mandelbrot fractal" on the pcanvas by examining behaviors of the orbits. Figure 1.3 shows a couple of examples.
Figure 1.3. Mandelbrot Fractals Painted on pCanvases
"Untitled"

 "Jellyfish Couple"

Now, suppose p is a fixed constant in the dynamical system while z_{0} is a variable. Then by the symmetric argument, we can talk about a zcanvas comprising pixels whose representative complex numbers are interpreted as values of z_{0}. Each pixel, er complex number z_{0}, on the zcanvas gives rise to a unique sequence z_{n} (with the fixed value of p) called the orbit of z_{0}, and as we shall see in § 6, we plot a fractal called a "Julia fractal" on the zcanvas by picking on certain behaviors of these orbits.
Figure 1.4. Julia Fractal Painted on zCanvas
There is a special subset of the Julia fractals, consisting of fractals generated in essence by socalled "Newton's rootfinding algorithm." The fractals in the subset are called "Newton fractals," which we will discuss in § 7.
Figure 1.5. Newton Fractal Painted on zCanvas


Back to the
Top


§ 2. The Divergence Scheme
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 gets larger. Our goal of § 2 is to introduce a fractal plotting technique associated with the notion of divergence of orbits of p using the aforementioned Mandelbrot equation (1.1) as an example.
Define a function f_{p} of a complex variable z involving a complex parameter p by setting f_{p}(z) = z^{2} + p and write the Mandelbrot equation 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  § 4, we set
(2.2) z_{0} = 0,
which is the critical point, so all orbits of p given by (2.1) are critical orbits, the orbits having the critical point as their initial values z_{0}.
Regard the complex plane as the set of parameters p of (2.1). Since a whole lot of the (critical) orbits of p in (2.1) appear to diverge to ∞, we are interested in the region in the complex plane comprising the parameters whose orbits do not diverge to ∞. One of the advantages in using the critical orbits (more of which will be seen in § 6) is that it allows us to prove the following theorem that turned out to be extremely useful.
The Divergence Criterion:
z_{m} > 2 for some m if and only if the orbit z_{n} of p diverges to ∞.
There are a couple of immediate corollaries to the theorem: First, the divergence criterion remains true if we replace 2 by any real number θ ≥ 2. Second, if p > 2, then z_{1} = p > 2 by (2.1), and consequently, the orbit z_{n} of p diverges to ∞; hence the region in the complex plane comprising the parameters whose orbits do not diverge to ∞ lies in the circle of radius 2 about the origin.
Suppose R is the square canvas bounded by xmin = 2, xmax = 2, ymin = 2 and ymax = 2 containing the circle of radius 2 and comprising 2,000 × 2,000 = 4,000,000 pixels. Regard R as a pcanvas so each pixel (i. j) in the canvas is identified with a unique representative parameter p belonging to it.
Let M = 1000 and θ = 2 and for every pixel (i. j), er parameter p, in the pcanvas R, iterate (2.1) at most M times and compute the nonnegative integer
(2.3) d(i, j) = the least (or first) iteration index m < M such that z_{m} > θ,
if such m exists, and d(i, j) = M, otherwise. Here, M is the maximum number of iterations which is intended for thwarting the computer to get trapped in an infinite loop. Also, in actual computation, we use z_{m}^{2} > θ^{2} instead of z_{m} > θ to avoid the "hidden" square root operation. Thus, each d(i, j) represents, unless d(i, j) = M, the smallest number of iterations m or "time" it takes for the orbit z_{n} of z_{0} to escape from the circle of radius θ before taking a long journey toward ∞.
The formula (2.3) generates an array d(i, j) of the least iteration indices, or an iteration array for short, over the canvas R, which can be converted to a visual image on the canvas. For example, if we paint the whole canvas white initially and then color the pixels (i, j) red if d(i, j) < M and is even, and color it black if d(i, j) < M and is odd, we get Figure 0.1. The image in Figure 2.1 on the left shows a cropped and resized version of the figure. The image on the right is similarly generated except that it uses θ = 10 instead of θ = 2 in (2.3).
Figure 2.1. The Mandelbrot Set
z_{m} > θ, θ = 2

 z_{m} > θ, θ = 10

The Mandelbrot Set: By definition, the famous Mandelbrot set is the set of all complex parameters p in the complex plane whose critical orbits do not diverge to ∞, or equivalently, the set of all parameters p whose critical orbits stay within the circle of radius 2 forever. The "snowman" of Figure 2.1 on the left, which is left untouched by the above redblack painting scheme and retains the initial white canvas color, depicts an approximation of the Mandelbrot set on the pcanvas given by replacing "forever" by "up until M = 1000 ".
The Divergence Scheme: We call the aforementioned process of computing an iteration array based on (2.3) the divergence scheme. In § 4, we'll see what we call the "convergence schemes" to deal with the orbits that "converge" to various "cycles."
Threshold Values: We call θ that appears in (2.3) a threshold value of the divergence scheme. Figure 2.1 shows that when θ is greater, the redblack pattern is more detailed and it is likely to produce better local images as well as a better approximation for the Mandelbrot set. On the other hand, it is intuitively clear that the computation is more timeconsuming when θ is greater, and it explains why Mandelbrot used the smallest threshold value available, namely θ = 2, when the computers were incomparably slower than today's counterparts.
There is another merit in using a larger threshold value. When we use a different dynamical system or a noncritical point for the initial value to explore wider fractal plotting possibilities, we often don't know what threshold value to use in the divergence scheme and a large threshold value is likely to provide a solution for it. But be careful: increasing the size of a threshold value quickly reaches a point of diminishing returns like increasing the maximum number of iterations for the accuracy of the approximation.
The Most Important Fact: Figure 2.1 shows that the redblack pattern in either image generated by the divergence scheme keeps getting more intricate without limit as we keep zooming in on an area nearer the boundary of the Mandelbrot set. It provides us with boundless possibilities of magnifying microscopic areas near the boundary and capturing fascinating images.
Example: If we call the images of Figure 2.1 global images of the Mandelbrot set, the image in Figure 2.2 shown below is a local image given by zooming in on a small neighborhood of a point that is just outside of the Mandelbrot set but very near its boundary.
It is a cropped and resized image from a figure on a large pcanvas with 6,400 × 3,200 pixels centered at the complex number (0.2820607, 0.011014375) with the horizontal radius 0.0000011. The part of the Mandelbrot set shown in the image is painted black and the exterior of the Mandelbrot set in multiple colors using the techniques described in Fractal Coloring. Note that the figure shows a distinct selfsimilar fractal structure as well as 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
We used the large canvas and M = 100,000 as the maximum number of iterations to generate the above image so as to get a possibility for zooming in further and capturing an additional image like Figure 0.6 as well as a possibility for making a high resolution printout of the image as an option.
.
§ 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. Although it may not sound obvious unless we know something about fractal dimensions, the following celebrated theorem guarantees that no figures on the plane are more complex 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.
But we wonder. Does it apply to, say, the image in Figure 2.2? Its colorful and intricate patterns surely look impressive, but does it have the maximum complexity implied by the theorem? We know that the boundary of the Mandelbrot set is somehow responsible for the striking patterns but it is colored black and mostly invisible in Figure 2.2.
It turned out that we can actually see what the boundary of the Mandelbrot set is like by lighting up its thin filaments:
Figure 3.1. The Boundary of the Mandelbrot Set in Figure 2.2
It shows that the Mandelbrot set in the rectangular area is vividly selfsimilar and composed of selfsimilar replicas of a large part of the image. Through the selfsimilarity, we can see it contains infinitely many replicas of the "snowman" painted black, although most of them are too small to be visible. The luminous boundary of the Mandelbrot set appears to be an elaborate network connecting these snowmen, and it gets so dense near the snowmen, it fills infinitely many areas of the plane like a spacefilling curve. The observation provides us with an intuitive idea as to why the "fractal dimension" of the boundary of the Mandelbrot set is equal to the topological dimension of the plane and why it is equated with the complexity of the Mandelbrot set.
The Mandelbrot set with its complex boundary shown above turns the complement (outside) of the Mandelbrot set into an extremely convoluted maze. Suppose we are shrunk to a pixel size and trapped in the maze painted in darker colors. Can we get out of the maze? The aforementioned selfsimilarity seems to show that the answer is affirmative even though the walls seem to be, quite amazingly, connected as one piece as well. The following fundamental theorem about the Mandelbrot set answers the question definitively. Here, a set is "simply connected" if it has no loop to trap anybody in it.
The DouadyHubbard Theorem: The Mandelbrot set is connected and simply connected.
Shown below on the left is another fractal generated by the Mandelbrot equation and the divergence scheme, where the Mandelbrot set is colored black and invisible. When its boundary's extremely thin "tungsten filaments" are turned on, we get the "nighttime" image on the right, vividly showing the presence of the Mandelbrot set invisible in the "daytime" image on the left.
Figure 3.2. Fractal in Daytime and Nighttime
Figure 3.2 again shows that fractals given by the Mandelbrot equation we normally see are daytime images like the image on the left and painted on the complement of the Mandelbrot set. The nighttime images certainly make it easier for us to visualize Shishikura's theorem and the DouadyHubbard theorem.
Plotting the complex boundary of the Mandelbrot set with reasonable accuracy may demand weeks of computing time even with a fast modern computer. Figure 3.2 shown below is a resized image from a fractal on the pcanvas with 4,000 × 4,000 pixels centered at the point
p = (0.25000316374967, 0.00000000895902)
with a microscopically small radius ≈ 0.0000000000001 = 10^{13}. We note that p is very near the cusp (0.25, 0) of the cardioid body of the snowman in Figure 2.1.
Figure 3.3. A Snowman Under the Microscope
M = 1,500,000

 M = 500,000

For the image on the left, we used whopping 1,500,000 iterations of the Mandelbrot equation for each black pixel because of ∞ involved in the definition of the Mandelbrot set. If we use M = 500,000 (still a large number) instead, the outline of the snowman becomes blurry as shown in the picture on the right. Fortunately, computers are inexpensive nowadays and we can afford a second or third computer to do tedious jobs. Programming carefully so as to minimize computing time is not as important as it used to be.
Shown below is the boundary of the Mandelbrot set in Figure 3.3, where it is invisible. It shows what was computed by the billions of extra iterations of the Mandelbrot equation in order to sharpen the blurry part near the boundary of the snowman.
Figure 3.4. The Boundary of the Mandelbrot set in Figure 3.3
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.2. 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 Convergence Schemes
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 "snowman" shape in Figure 2.1 by heart: The main body is in the form of a heartshaped "cardioid" with a bunch of circular disks attached, and to each of the disks another bunch of disks are attached. The pattern repeats as if the cardioid has children, grandchildren, great grandchildren and so on and so forth, but beyond that nobody knows exactly what's happening. Here, the cardioid means the bounding curve together with its interior.
As we have seen, the divergence scheme paints the interiors of the cardioid and all of the disks in a single color like white as it is incapable of distinguishing them. Our current goal is to find their distinguished mathematical properties and paint them in various colors, like Figure 4.1 below, by developing a scheme that is different from the divergence scheme.
Figure 4.1. The Interior and Boundary of the Mandelbrot Set
A sequence c_{n} of complex numbers is called a cycle if there is a positive integer k satisfying c_{n} = c_{n+k} for any index n. The smallest such integer k is called the period of the cycle, and a cycle with period k is called a kcycle for short. For example, the sequence
1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, · · ·
is a 3cycle but not a 6cycle or a 9cycle. The sequence 0, 0, 0, 0, · · · is a 1cycle, which we identify with the constant 0.
A sequence z_{n} is said to converge to a kcycle c_{n} if it gets closer and closer to c_{n} or behaves more and more like c_{n} as n gets bigger. If a sequence converges to a kcycle, we also say that the sequence is attracted to the kcycle.
Example 1: The sequences 1/2, 1/3, 1/4, 1/5, 1/6, · · · and 1/2, 2/3, 3/4, 4/5, 5/6, · · · converge to the constants 0 and 1, or equivalently, to the 1cycles 0, 0, 0, 0, · · · and 1, 1, 1, 1, · · ·, respectively. Therefore, the sequence
1/2, 1/2, 1/3, 2/3, 1/4, 3/4, 1/5, 4/5, 1/6, 5/6, · · ·, 1/1000, 999/1000, 1/1001, 1000/1001, · · ·
converges to the 2cycle 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, · · ·, 0, 1, 0, 1, · · ·.
Suppose ε is a small real number and z_{n} an orbit of parameter p as in § 2. We note that if z_{n} behaves more and more like a kcycle c_{n} as n gets bigger, then because c_{n+k}  c_{n} = 0, we surely have z_{n+k}  z_{n} < ε for large n. This gives us:
Proposition: If the orbit z_{n} converges to a kcycle then z_{n+k}  z_{n} < ε for sufficiently large n.
The proposition provides us with our convergence scheme with period index k, which is given by replacing z_{m} > θ of the divergence scheme by the inequality z_{m+k}  z_{m} < ε. Thus, the convergence scheme generates the iteration array
(4.1) d(i, j) = the least (or first) iteration index m < M such that z_{m+k}  z_{m} < ε
over the canvas R comprising pixels (i, j), which in turn can be converted to a visual image by fractal coloring. In (4.1), z_{n} is the orbit of p which is the parameter representing the pixel (i, j), M is a maximum number of iterations (just like in the divergence scheme) and ε is a small positive real number such as 10^{6} or min{Δx, Δy}; see (1.2). There are no set formulas for M and ε so we find them by educated guesses through our experience and computer experiments.
The convergence scheme with period index k is intended not only to compute the iteration array d(i, j) but also to detect every parameter whose orbit converges to a kcycle. We will talk about "a parameter whose orbit converges to a kcycle" frequently, so let's call it a parameter of period k for short. In particular, a parameter of period 1 is the same as a parameter whose orbit converges to a single point (1cycle).
Figure 4.2. The Mandelbrot Set by the Divergence and Convergence Schemes
Example 2: The images in Figure 4.2 are given by applying the divergence scheme and the convergence schemes with the period indices k = 1, 2, · · ·, 96 to the Mandelbrot equation. The convergence scheme with period index k = 1 detects that the interior of the cardioid of the Mandelbrot set comprises parameters of period 1 and paints it in various shades of red, and the convergence scheme with period index k = 2 detects that the interior of the largest disk comprises parameters of period 2 and paints it in orangish colors. The repetitive use of the convergence schemes up to k = 96 with various colors generates the image on the right.
Figure 4.3. Closeups of the Mandelbrot set by the Convergence Schemes
Note that the small replicas of the snowman that look like isolated islands are actually connected by the boundary of the Mandelbrot set; see Figure 4.1.
Periodicity Diagram: If the interior of, say, a disk in the Mandelbrot set comprises a parameter of period k, let's say the disk has period k. For example, we have seen that the largest disk has period 2. If we label the cardioid and disks of the Mandelbrot set with its periods 1, 2, etc., instead of coloring red and orange, etc., we get Wikipedia's Periodicity Diagram instead of the colorful image shown above. The periods in the diagram have interesting numerical patterns that are easy to recognize and will play an important role when we plot colorful "Julia fractals" in § 6. The numerical patterns are yet another amazing property of the Mandelbrot set.
Weakness of the Convergence Schemes: The convergence scheme is based on the aforementioned proposition whose converse is unfortunately not necessarily true. For example, the convergence scheme with k = 1 may be fooled by an orbit that diverges to ∞ very slowly and treat it like a convergent orbit. Also, the convergence scheme with, say, period index 6 cannot distinguish parameters of periods 1, 2, 3 and 6 and we can locate parameters of period 6 only after identifying all parameters of periods 1, 2, 3. For most of the practical problems, however, we can overcome these shortcomings by carefully programming a computer to do the right thing.
Brute Force Algorithm: Because of the aforementioned weakness, we can compute the iteration array for an image like Figure 4.2 on the right only by the timeconsuming brute force algorithm that goes like this: For each pixel (er parameter) on the canvas, check if its orbit diverges using the divergence scheme; if it does, compute d(i, j), but if not, check if it has period 1 using the convergence scheme with period index 1; if it does, compute d(i, j), but if not, check if it has period 2, etc. etc. Repeat the process until it comes to the last period index k like k = 96 for the image of Figure 4.2 on the right. The algorithm is simple and straightforward, but it is often called an "exponential time algorithm" as its computing time grows almost exponentially with k.
The Eyeball Effect: Weaknesses and shortcomings are not all that bad and sometimes they can give us unexpected artistic effects. Even a programming bug can give us beautiful output sometimes. The picture below on the left is essentially the same as Figure 3.3 and is given by the divergence scheme alone while the one on the right is painted by the divergence scheme (with different colors) and the convergence scheme with period index k = 1. The "eyeballs" painted by the convergence scheme are caused by its weakness of confusing some of the slowly divergent orbits as convergent. The images show which parameters are affected. Figure 6.8 illustrates the "eyeballs" more vividly.
Figure 4.4. The Eyeball Effect (Right) Given by the Convergence Scheme
§ 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, a 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 elementary function containing a parameter p. The initial value for the orbits is often a critical point of f_{p} as in the case of the Mandelbrot set but it is not a requirement. There are functions f_{p} without convenient critical points, and even if f_{p} has one, a noncritical point often produces interesting fractals from f_{p}. We always have the advantage of having a fast computer that allows us to satisfy our curiosity by engaging in a computer experiment.
Figure 5.1. "Snowman" Shot by Arrow
The comical image shown above is a Mandelbrot fractal given by
(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 shot the arrow); see the inset for an enlarged image. Here, the picture is actually given by rotating the original output 90^{o} counterclockwise to better fit in the webpage. In this article, we use rotations and horizontal/vertical reflections of fractal images freely for artistic effects.
The image shown below is a Mandelbrot fractal given by zooming in on the boundary of the tiny snowman. We often use it as a night sky of 3D landscapes like in Figure 5.2 and Figure 5.3 shown below. § 8 explains how to plot such 3D images.
Figure 5.2. "Mandelbrot Island"

 Figure 5.3. "Desert Rocks"

The Logistic Equation: In 1838 Pierre Verhulst introduced a differential equation called the "logistic equation," which became a widely used mathematical model for population dynamics. If we replace the derivative in the equation by its approximating difference quotient and do some algebra, we get the following formula that is more suitable for computer applications:
(5.3)
z_{n+1 }= f_{p}(z_{n}) = p(1  z_{n}) z_{n }.
It is again called the logistic equation (or logistic map) and is equally applicable in the population dynamics when the variables are restricted to real numbers.
In 1974, while conducting a computer simulation of certain population change, biologist Robert May discovered "very complicated orbits" of (5.3), which led us to the concept of chaos. In 1993, a "chaotician" appeared in Steven Spielberg's hit movie, "Jurassic Park," tacitly suggesting a possibility of chaos in the controlled dinosaur populations.
Figure 5.4. "Bifurcation"
So, it is natural that we plot Mandelbrot fractals of the dynamical system (5.3) by expanding its variables to complex numbers. Figure 1.1 and Figure 5.4 shown above provide examples. The horizontal line of Figure 5.4 is a part of the real interval [α, 4] with α ≈ 3.57 on which Robert May discovered "chaotic" orbits. In this example, a noncritical point z_{0} = 0.1 is used for a slight deforming effect.
By the way, people familiar with multivariable calculus might find a fun project in mapping an image like Figure 5.4 on surfaces like the ones shown below. We can see more examples of this sort in Gallery 3D.
The next two images are Mandelbrot fractals given by the cubic equation motivated by the Logistic Equation:
(5.4)
z_{n+1 }= f_{p}(z_{n}) = p(1 + z_{n})(1  z_{n}) z_{n }.

 Figure 5.5. "Rising Dragon"

§ 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 comprising parameters p. A fractal is called a Julia fractal if it is given by a dynamical system of the same form with a fixed parameter p and painted on a zcanvas comprising initial values z_{0} instead. It is named after Gaston Julia, who was one of the early pioneers of fractals generated by dynamical systems of complex numbers.
Thus, if we use the divergence scheme to generate a Julia fractal, then for each pixel (i, j) in the canvas R, we use the complex number representing the pixel as an initial value z_{0} and iterate (5.1) with the fixed parameter p to compute the orbit z_{n} of z_{0} until we get
(6.1) d(i, j) = the least (or first) iteration index m < M such that z_{m} > θ,
if such m exists, and d(i, j) = M, otherwise. As we have seen earlier, M is a maximum number of iterations and θ is a threshold value. If we use the convergence scheme with period index k, then instead of (6.1), we have
(6.2) d(i, j) = the least (or first) iteration index m < M such that z_{m+k}  z_{m} < ε,
if such m exists, and d(i, j) = M, otherwise. Here, ε is a small real number such as 10^{8}. d(i, j) collectively forms an array over the canvas R called an iteration array, which in turn can be converted to a visual fractal by Fractal Coloring.
Let's follow various concepts related to Julia fractals with the familiar Mandelbrot equation (again)
(6.3) z_{n+1 }= f_{p}(z_{n}) = z_{n}^{2} + p.
Example 1: Figure 6.1 shows two nearly equal parameters that are fixed and two Julia fractals given by (6.3) and the respective parameters. The first parameter belongs to the interior of the Mandelbrot set, or more precisely, the interior of a disk of period 11 and the second parameter to the outside of the Mandelbrot set. We call the Julia fractals "Hydra with Eleven Heads" and "Hydra's Ash," respectively. The green/blue "Hydra" on the left is given by the convergence scheme with period index k = 11 and its background by the divergence scheme. "Hydra's ash" is given by the divergence scheme alone.
It is another fascinating fact about the Mandelbrot set that the period of a disk where the parameter p resides always has a geometric implication in the Julia fractal such as the number of heads in the figure, but why it is so is not completely understood. If the period is a large composite number rather than a prime and the disk is not directly attached to the "big caridoid" of the Mandelbrot set, the outcome gets really interesting.
Figure 6.1. "Hydra with Eleven Heads" and "Hydra's Ash"
Julia Sets: Recall that the Mandelbrot set is by definition the set of all parameters p in the complex plane whose critical orbits do not diverge to ∞. As we have seen, we can always view the complex plane as the set of parameters p or the set of initial values z_{0} for the orbits z_{n} with a fixed parameter p.
If p is a fixed parameter, then by the filledin Julia set of p, we mean the set of all initial values z_{0} in the complex plane whose orbits z_{n} (with the fixed p) do not diverge to ∞. Thus, the definitions of a filledin Julia set and the Mandelbrot set are similar, but unlike the Mandelbrot set, there are infinitely many filledin Julia sets. Like the Mandelbrot set, we can approximate many of the filledin Julia sets by using the divergence and/or convergence schemes and plotting them on a canvas with finitely many pixels. For example, "Hydra with Eleven Heads" of Figure 6.1 without the reddish background is a computer approximation of a filledin Jula set.
By the Julia set of p, we mean the boundary of the filledin Julia set. As the names suggest, the filledin Julia set intruoduced earlier plays a supporting role for the Julia set in mathematics, and it is because a Julia set is where "chaos" occurs as well as it is a fractal whose fractal dimension varies and may equal to that of the boundary of the Mandelbrot set. Julia sets are among the most important in fractal geometry, but in fractal art, it is usually filledin Julia sets we look at and appreciate. It is often difficulit to plot a Julia set alone without its filledin inside, but if we manage to plot it, it can be stunning because of its complexity; e.g., see Figure 6.4.
The concept of Julia set naturally extends to a more general dynamical system (5.1), but a lot of things about it are still in mystery and belong to experimental mathematics by the use of computers. Here's one fascinating and useful fact however. The following theorem was established about a hundred years ago before the computer era by Gaston Julia and
Pierre Fatou independently and explains why the critical orbits are important.
The FatouJulia Theorem: Consider a dynamical system of the form
(6.4)
z_{n+1 }= f_{p}(z_{n}) = c_{m} z_{n}^{m} + c_{m1} z_{n}^{m1} + · · · + c_{2} z_{n}^{2} + c_{1} z_{n} + p ,
where m is an integer ≥ 2 and c_{m}, c_{m1}, · · ·, c_{2}, c_{1} are complex constants. Then the Julia set of (6.4) of p is connected if and only if every critical orbit of p stays within a finite bound.
Note that if m > 2 in (6.4), p may have multiple critical orbits as f_{p} may have multiple critical points. In case of the Mandelbrot equation (6.3), each p has a single unique critical orbit; hence the FatouJulia theorem is stated beautifully as:
The Julia set of p is connected if and only if p belongs to the Mandelbrot set.
For the Mandelbrot equation, it is also known that if the Julia set of p is disconnected then it must be a set of "totally disconnected" points called a Cantor dust (named after Georg Cantor, the pioneer of set theory) and cannot be the disjoint union of, say, three connected pieces. For example, because of our choice of the parameters in Figure 6.1, the Julia set (hence the filledin Julia set) in "Hydra with Eleven Heads" is connected (as one piece) while the Julia set in "Hydra's Ash" is a Cantor dust. The boundaries of the "Gold Lions" in Figure 6.2 shown below are also examples of connected Julia sets. Note that the connectedness is not entirely obvious if we rely only on our eyes.
Figure 6.2. "Gold Lions" of Periods 14 and 42 by the Mandelbrot Equation
Plotting Colorful FilledIn Julia Sets of the Mandelbrot Equation: In the following, we plot a few filledin Julia set using parameters p that are inside the Mandelbrot set so the Julia sets are connected according to the JuliaFatou theorem. Since the Mandelbrot set lies in the circle of radius 2 about the origin in the complex plane, suppose p ≤ 2. Then it is not particularly difficult to prove that any θ ≥ 2 works as a threshold value for our divergence scheme (6.1). The following shows our method to draw a connected Julia set of p:
Step 1: Use Wikipedia's Periodicity Diagram to choose a disk D of period k of our interest.
Step 2: Identify the disk D in our own computer plot of the Mandelbrot set and choose a pixel (parameter) p in the interior of D.
Step 3: Convert the pixel coordinates p = (i, j) to the Cartesian coordinates p = (x, y) using (1.3).
Step 4: Choose a zcanvas centered at z_{0} = 0.
Step 5: Use the divergence scheme and the convergence scheme with period index k to plot the filledin Julia set of p on the zcanvas of Step 4. The divergence scheme generates the Julia set together with its exterior and the convergence scheme the interior of the Julia set.
Note: We have to have the computer program for plotting the Mandelbrot set of Step 2 so we have the numbers necessary to do the conversion in Step 3. In Step 5, the convergence scheme with period index k plots the filledin Julia set of p whose boundary is the Julia set. If we want to plot the Julia set that is a Cantor dust, choose p in the exterior of the Mandelbrot set which is very near the disk of Step 2, follow Steps 3 and 4 and use the divergence scheme alone in Step 5.
Example 2: For the "Gold Lions" of Figure 6.2 shown above, note that there is a disk, call it D_{14}, of period 14 attached to the cardioid near its cusp and a disk D_{42} of period 42 = 14 × 3 attached to the disk of period 14 in the Periodicity Diagram.
For the
"Gold Lion" on the right, choose p = (0.296555, 0.020525) in the interior of D_{42} using Steps 2 and 3 and plot a Julia fractal using Steps 4 and 5. The "Gold Lion" is the filledin Julia set of p, which is connected as its boundary, the Julia set of p, is connected according to the FatouJulia Theorem, and it clearly shows a geometric pattern associated with the numbers 14 and 3. The left "Gold Lion" is given by p = (0.296498, 0.020525) chosen from the interior of D_{14}. Because the two parameters are almost equal, the two figures are similar, but the left "Gold Lion" does not show the number 3 as clearly as the other "Gold Lion." Recall that we use rotations and horizontal/vertical flippings of fractals freely for artistic effects.
Example 3: The "Red Lion" of period 68 = 17 × 4 of Figure 6.3 is similarly generated but it has a shape that is much more intricate than the "Gold Lions" shown earlier. That's because the parameter of the Julia set is chosen from very near the boundary of the disk of period 68. It is known that the fractal dimension of the Julia set of a parameter p gets closer to 2 if p gets closer to the boundary of the Mandelbrot set, and as we saw in § 3, it means that the complexity of the Julia set approaches the maximum complexity available for an image depcited on a 2D plane. In any event, the golden rule for getting an intricate fractal we saw in § 2 still holds for the Julia sets: Choose a parameter that is as close to the boundary of the Mandelbrot set as possible.
Figure 6.3. "Red Lion" of Period 68 = 17 × 4
The "Red Lion" painted in various shades of red is a filledin Julia set but its boundary, which is the Julia set, is colored black and nearly invisible in Figure 6.3. If we turn on its extremely thin "tungsten filaments" however, it suddenly comes alive as the following image shows:
Figure 6.4. The Julia Set (Boundary) of "Red Lion"
Julia Sets by Other Dynamical Systems:
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. There p = (0.185, 0.00007666), and the gold "Twin Dragons" is painted on a zcanvas centered at the origin by the convergence scheme with period index k = 2. If we alter the value of p and find its period, we get a variety of "Twin Dragons" such as:
Figure 6.5. "Twin Dragons" by the Dynamical System (5.2)
And here is another with a shining Julia set:
Figure 6.6. "Obese Seahorse Family" by the Dynamical System (5.4)
Similarities of Mandelbrot and Julia Fractals: Although the precise reason is unknown, we note that in many cases Mandelbrot and Julia fractals from the same dynamical system have similar appearances "locally" if they are generated by nearly equal parameters. For example, the image below on the left is a Mandelbrot fractal given by the logistic equation (5.3) on a microscopic neighborhood of p = (3.000008596, 0.076065598) and the one on the right is a Julia fractal from the same equation with the fixed parameter q = (3.001, 0.075975) ≈ p.
"Fractal Elephants" by the Logistic Equation
The image on the left in Figure 6.7 below is a Mandelbrot fractal given by the Mandelbrot equation (6.3) on a pcanvas centered at p = (0.25000316374967, 0.00000000895902) that belongs to the body of the "snowman." It is essentially the same as the images in Figure 3.3 whose complexity was discussed in § 3.
The image on the right is a Julia fractal given by the fixed parameter p on a zcanvas centered at the origin. Both images contain infinitely many "cuttlefish" with their eyeballs given by the eyeball effect of the convergence scheme. The divergence scheme painted all other areas except for the black "holes." Just like the Mandelbrot set discussed in § 3, the filledin Julia set in the Julia fractal on the right comprises infinitely many black holes, which is, according to the FatouJulia Theorem, connected by an invisible network and it shapes the colorful and intricate Julia fractal. The example shows that a Julia set can be as complex as the Mandelbrot set "locally."
Figure 6.7. "Partying Cuttlefish" by the Mandelbrot Equation
Figure 6.8. Closeup of "Partying Cuttlefish"
§ 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 the parameter 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 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 the NewtonRaphson RootFinding Algorithm, aka
Newton's Method. Hence, each orbit of (7.1) converges to a root of g quickly more often than doing something else, and it allows us to plot most of the Newton fractals by the convergence scheme (with period index k = 1) alone with a relatively small maximum number of iterations like M = 500.
Furthermore, if we know all the roots of g prior to the fractal plotting, we can modify the convergence scheme fairly easily so as to add more colors to Newton fractals of g. Because a Newton fractal is a Julia fractal, a "canvas" and an "orbit" always mean a zcanvas and an orbit of z_{0}, respectively, in this section.
Example 1 (Roots of Unity): Probably the simplest Newton fractals to plot are given by a polynomial of the form
g(z) = z^{ d}  1
as its roots are readily available by hand calculations or Googling "roots of unity."
Figure 7.1. Newton Fractals of g(z) = z^{ 5}  1
"Crab Queue"
The leftmost image of Figure 7.1 is a Newton fractal for g(z) = z^{ 5}  1 painted on a square canvas centered at the origin with radius 1.1. It uses five essentially different colors, sky blue, purple, red, amber, and blue, associated with the five roots of g. The sky blue region, e.g., comprises the initial values z_{0} in the canvas whose orbits converge to the root r = (1, 0) and is called the basin of attraction of Newton's method for the root.
Thus, there are five basins of attraction in the leftmost fractal and they are divided by the basin boundary. The basin boundary is precisely the Julia set of the Newton fractal, and that is where Newton's rootfinding algorithm behaves in a "chaotic" fashion and fails to provide a root. As we can see clearly, a Newton fractal is totally uninteresting unless it contains a part of the Julia set—just like the Julia fractals discussed in the preceding section. It is known that the Julia set is Cantor dust.
The second image of Figure 7.1 is a variation of the first and the third is given by zooming in on one of the "bands" in the second image.
Example 2 (Cyclotomic Polynomials): Another interesting example with known roots is a
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.
Figure 7.2. Cyclotomic Polynomials with Eight Roots
Finding the Roots: As we have seen, a Newton fractal of a polynomial with colorful basins of attraction requires its roots to be known. So, if we don't know the roots, how can we find them? A natural choice seems to be the use of Newton's method, but unfortunately, its chaotic nature makes it difficult to program a computer and consistently find the roots. Another way is to use Müller's Method instead. Although Müller's method lacks the impressive simplicity and speed of Newton's method, it generally works well and automatically finds all roots of the polynomial. All Newton fractals shown in Gallery 2D use Müller's method even when the polynomials have known roots.
Figure 7.3 is a Newton fractal of a fifth degree polynomial whose roots are given by Müller's method. Just for fun, we painted it on a plane, a sphere and a torus.
Figure 7.3. Newton Fractal on Plane, Sphere and Torus
Figure 7.4 shows a Newton fractal of a 12^{th} degree polynomial painted on a plane and a sphere. The second image which is on a sphere is intended to give a 3D appearance. All of the twelve roots are again found by Müller's method quickly.
Figure 7.4. "Dragonfly" on Plane and Sphere
Figure 7.5 illustrates two Newton fractals of a fifth degree polynomial. We painted them by highlighting different parts of essentially the same image. Note that we can find numerous replicas of the "fish" in the "crab" and vice versa. Note also that each of the images shows white Cantor dust which is a part of the Julia set (or basin boundary) of the Newton fractal.
Figure 7.5. "Newton Fish" and "Newton Crab"
