High-Order Mandelbrot and Julia Sets

High-Order Mandelbrot and Julia Sets
Written by Christopher Thomas

Normal Mandelbrot and Julia Sets

Most people know about the Mandelbrot set and Julia sets nowadays. They were among the first fractals whose pictures touched public awareness; their intricate complexity decorates many walls and computer monitors.

The Mandelbrot set and Julia sets are both produced by iterating the following equation using complex numbers (numbers with real and imaginary components, representing points on the Argand plane):

Z <- Z² + C

Iterating this equation causes a point to hop around on the plane. The point will either end up staying near the origin, or will fly off to infinite distance. From a derivation which will be shown later, we know that if the point ever gets outside a radius-2 circle about the origin, it is forced to fly out to infinity.

A point is inside the set, then, if it never goes outside the radius-2 circle. We can produce pretty pictures by tracking points that go outside the circle, and colour-coding them based on the number of steps it took to leave the circle. Points that stay in are assumed to be part of the set (some won't be, but it would take an infinite amount of time to tell for sure using this method).

The same equation is used to produce both Mandelbrot and Julia sets. To draw a Julia set, we select a fixed value of C and use the location on the plane as the starting value of Z. To draw the Mandelbrot set, we use the location on the plane as the value of C, and set Z to zero to start with (or set it to C, which is what it ends up being after one step, but if you think of it as zero it makes the next paragraph make more sense).

When you think about it, there are actually four variables to be set when initializing the equation - the real and imaginary components of both Z and C. This means that you can think of the Mandelbrot and Julia sets as all being cross-sections of a four-dimensional object. The Julia sets are the set of all cross-sections parallel to one pair of axes. The Mandelbrot set is one specific cross-section parallel to the other pair of axes (the plane defined by "Z = 0"). You can get strange Mandelbrot-like sets by setting the initial value of Z to a different value. These, along with the Mandelbrot set, are another complete set of cross-sections of the four-dimensional object.

Mandelbrot Set with Z₀ = 0.5+0.5i

Three-dimensional cross-sections of this object look very neat, but are beyond the scope of this discussion.

High-Order Mandelbrot and Julia Sets

We now consider Mandelbrot and Julia set analogues defined by a similar equation with an exponent other than 2:

Z <- Zⁿ + C

The case of integer n is the easiest to implement. The iterated equation is easy to implement, and the only other algorithmic change is to change the radius of the limiting circle. For reasons discussed in the next section, the limiting radius is defined as:

r_limit = 2^{1 / (n - 1)}

That is, the (n-1)th root of 2.

Some pictures of higher-order Mandelbrot and Julia sets are below:

n=3 Mandelbrot Set	n=3 Julia Set with C = -0.5+0.5i	n=3 Julia Set with C = 0+1i
n=7 Mandelbrot Set	n=7 Julia Set with C = -0.621+0i	n=7 Julia Set with C = 0.1+0.95i

Some patterns emerge. Julia sets with order n have n-fold rotational symmetry, while Mandelbrot sets with order n have (n-1)-fold rotational symmetry. The reasons for this lie in the symmetries of exponentiation in the complex plane (described in the next section). The full proof of Mandelbrot and Julia set symmetries, however, is beyond the scope of this discussion.

Non-Integer Values of n

Non-integer exponentiation of real numbers is defined in terms of the exponential function:

a^b = e^{b ln a}

The same formula can be used with complex values of a and b, when we express a as an exponential using the Euler relation:

e^ki = cos(k) + i sin(k)

This gives the following relation in the general case:

e^{a + bi} = e^a (cos(b) + i sin(b))

Which is to say, the value of a is the natural logarithm of the radius of the point in polar coordinates, and the value of b is the angle of the point (in radians) in polar coordinates. By converting the location of a point on the complex plane into polar coordinates, we can therefore turn it into an exponential. This lets us take the logarithm and plug it into the exponentiation equation, and then we can apply the general version of Euler's relation again to find the exponential of the complex number.

The upshot of all of this is that, in polar coordinates, raising a to the power of b multiplies the angle by b, and raises the radius to the power of b:

c = a^b
c_r = a_r^b
c_angle = b x a_angle

This lets us evaluate the equation Z <- Zⁿ + C for non-integer (and non-real) n.

It turns out, however, that there are problems. The relation

c_angle = b x a_angle

has more than one solution for non-integer values of b. For rational values of b, it has a finite number of solutions; for irrational, infinite. In neither case is the correct interpretation of the result clear. Detailed exploration of the various interpretations, and of complex values of b, is beyond the scope of of this discussion.

Lastly, the polar coordinate interpretation of exponentiation explains how the limiting radius is derived. For an exponent of n, implementation of Z <- Zⁿ gives a maximum new radius of

r_new = r_oldⁿ

In the limiting case, the largest possible value of C applied in the best possible direction will bring the value of r_new back to r_old. In other words, the pull away from the origin is exactly balanced by a displacement back towards it. Given that the largest acceptable value of C is also the limiting radius, and subbing in r_limit for r_new and r_old, we get:

r_limit = r_limitⁿ - r_limit
2 x r_limit = r_limitⁿ
2 = r_limit^{(n - 1)}
2^{1 / (n - 1)} = r_limit

Source Code

The program used to generate the images in this document is available here.

Known problems are as follows:

Handling of non-integer n is incorrect.

Only one solution is picked, which means that sets with rational n just look like twisted versions of sets with integer n with discontinous slices, instead of shapes with highly increased symmetry.
The iteration equation implementation is slow and error-prone.

The iteration equation implementation converts to polar coordinates for the exponentiation, and then back to cartesian for the addition. This lets it handle almost any possible value of n, but with non-integer n failing to work properly, it's a waste of computing power. It also makes rounding errors the worst of both worlds, and introduces precision errors if the implementation of trig functions isn't perfectly accurate.

The pictures are still very pretty, though, and are correct for integer n.

Written by Christopher Thomas. Last updated 21 March 2003.