Gaston Julia

I always wanted to learn more about complex numbers. Not just the fact that the root of -1 was equal to i. I was particularly interested in figures with fractal structure. I wanted to understand what this meant and how to make such visualization. Having some time from work, I found a book on the theory of functions of a complex variable on the shelf and decided that it was time to start.

Iterated Function

Before writing an algorithm, we’ll have to learn a few basic notions. We’ll start with the concept of the iteration of function f. Let’s introduce the following notion of the n-th iteration of function f:

Iterated function

in which ○ is a function composition.

We’ll take the following identity map as the zero iteration:

Identity map

There are plenty of fixed-point theorems of certain types of maps. Let’s recall what a fixed point of the g map is. It’s a solution of the following equation:

Fixed-point theorem

For the iteration of the function, we also introduce the notion of a fixed point. It is a point from the domain of function f. Thus, the sequence of values of function f iterations converges to a fixed point of the f map:

Julia set. Fixed point.

in which y is a point sufficiently close to the x point.

The set of all y, satisfying the above condition is called the limit set of point x under the f function iteration.

As a rule, we’re interested in studying the sequence of functions formed by iterations (you can find a similar idea in Newton's method; the sought solution is the attracting fixed point of a built map and, therefore, can be found as the limit of a sequence of iterations).

Julia set. Sequence of functions.

as well as the sequence of values computed in some starting point

Julia set. Sequence of values.

It’s a bit problematic to visualize the iterative convergence of point y to x under the f function iteration for the complex world, as it is 4D. But for the real world, the picture is only 2D (the red trajectory diverges, while the blue one converges to a fixed point).

Julia set. Iterative convergence.

Thus, we have laid the first corner stone for understanding the Julia set. Let’s move on to the next one.

Complex Quadratic Polynomial

In general, we make Julia sets for the function of the following form:

Function to build Julia sets

in which p and q are complex polynomials.

To simplify the task, we will make the Julia set for a quadratic polynomial only. In general terms, a complex quadratic polynomial looks the following way:

Julia set. Complex quadratic polynomial.

There is also a concept of unitary (a = 1) centered (b = 0) form of a quadratic polynomial.

Julia set. One of the quadratic polynomial forms

Let’s introduce the following φ map, and find its inverse mapping:

Take a look at the following expression:

It is easy to see that this expression reminds a quadratic polynomial in the general form: p(z). Let’s find a c to convert the obtained expression to a quadratic polynomial of the general form:

Thus, we have shown that the above replacement of c allows to write a polynomial of general form like this:

Quadratic Polynomial Behavior under Iterations

Let’s move on to the third corner stone for understanding the Julia set. Let’s consider the second iteration of the complex quadratic polynomial:

Julia set. Complex quadratic polynomial.

It is not difficult to show that

Let’s formulate the conclusion: to study the behavior of a quadratic polynomial under iterations, it’s enough to study it in the unitary centered form, rather than in its general form. Actually, it’s really cool, as there are three coefficients in the general form, while there’s just one in the unitary centered form.

Chaotic and Normal Behavior

We will be able to give the definition of the Julia set pretty soon. We will consider intuitive definitions, rather than strictly mathematical ones. But first, let’s review the concept of stability. The solution of a differential equation is called stable if the solution behavior with conditions, close to the initial ones, is not much different from the behavior of the initial solution. It’s a no-brainer that the whole thing is about the phrase “is not much different”. Generally, there are three types of stability. For instance, Lyapunov stability is involved in the definition of the Julia set (this definition follows from Montel's theorem), but we will not go into the details of different stabilities, as the main thing is to get the idea. Let's try to illustrate this. To begin with, look at the intuitive perception of stability. After small changes, blue and yellow trajectories described by points behave always in the same way, while even the slightest change moves the red trajectory in a completely different direction.

Let’s consider the f(z) = z2 — 1 + 0.213i polynomial and build a trajectory for 100 iterations at

The above example shows that a small change in the initial condition resulted in a divergence of trajectories. In the area of the starting point, we will call such behavior chaotic. That is to say, behavior depends much on small changes in initial conditions. As for the behavior, during which stability is preserved, i.e. small changes in the initial conditions do not strongly affect behavior, we will call it normal, like in the following example. Look at the f(z) = z2 — 1 + 0.2i polynomial (we changed the c constant a bit) and build a trajectory for 100 iterations, when 

Fatou, Julia and Mandelbrot Sets

Attractor’s Basin of Attraction

So, we have learnt the meaning of the Julia set. Now we can get down to the algorithm of constructing such a set. But we should first introduce some definitions.

Let’s consider A(∞) in the context of our discussions for some f(z) = z2 + c polynomial.

Turns out, there is the following theorem. The A(∞) set (in the above context) is open, connected and unbounded. It is completely contained in the Fatou set. The Julia set coincides with the boundary of A(∞), which is a closed and bounded subset of all complex numbers.

Making the Julia Set

Reminding you that we will write an algorithm to make the Julia set of quadratic dynamics, not any function. Thus, we have learnt that the Julia set coincides with the boundary of A(∞), which is a closed and bounded subset of all complex numbers, i.e. the Julia set is also closed and bounded. Let’s denote a subset of complex numbers that is generated by iterations of function f and remains bounded. We will call it the filled Julia set:

Filled Julia set

The last theorem we need sounds the following way:

In addition to the algorithm, we can make a conclusion that the entire Julia set is inside a sphere of radius R with the centre at the origin of coordinates.

P.S. that’s exactly when the fans of the golden ratio, design, mysticism and occultism will notice that when f(z) = z2+ 1, threshold R is equal to the golden ratio о_0.

The Algorithm

  1. Select с to define the f(z) = z2 + с polynomial
  2. Compute R for the defined f(z) = z2 + с polynomial
  3. Select the maxIter parameter to denote the maximum iteration. Obviously, the higher it is, the better the accuracy, the slower the algorithm is.
  4. Generate an array of colors. The size is no more than maxIter. The change from a less bright color to a brighter one will denote how far the point is located from the Julia set.
  5. For each point, we determine whether it is a part of the filled Julia set, and also the number of iteration, during which the threshold has been exceeded.
  6. If |z| > R, we use the first color. Then, we use the color, at the iteration number of which the threshold has been exceeded.

The Procedure for Generating the Julia Set

/// <summary>
/// Generate bmp with Julia set
/// </summary>
/// <param name="c">parameter of square dynamics</param>
/// <param name="w">width of bmp</param>
/// <param name="h">height of bmp</param>
/// <param name="maxIter">max iterations of function</param>
/// <param name="xMin">window in complex plane</param>
/// <param name="yMin">window in complex plane</param>
/// <param name="xMax">window in complex plane</param>
/// <param name="yMax">window in complex plane</param>
/// <returns></returns>
static XBitmap PlotJuliaSet(ComplexNumber c, int w, int h, int maxIter,
    double xMin = Double.NaN, double yMin = Double.NaN, double xMax = Double.NaN, double yMax = Double.NaN)
    double r = CalculateR(c);
    if (Double.IsNaN(xMin) || Double.IsNaN(xMax) || Double.IsNaN(yMin) || Double.IsNaN(yMax))
        xMin = -r;
        yMin = -r;
        xMax = r;
        yMax = r;
    Logger.Instance.Log("R = " + r);
    double xStep = Math.Abs(xMax - xMin) / w;
    double yStep = Math.Abs(yMax - yMin) / h;
    XBitmap bmp = new XBitmap(w, h);
    IDictionary<int, IDictionary<int, int>> xyIdx = new Dictionary<int, IDictionary<int, int>>();
    int maxIdx = 0;
    for (int i = 0; i < w; i++)
        xyIdx.Add(i, new Dictionary<int, int>());
        for (int j = 0; j < h; j++)
            double x = xMin + i * xStep;
            double y = yMin + j * yStep;
            ComplexNumber z = new ComplexNumber(x, y);
            IList<ComplexNumber> zIter = SqPolyIteration(z, c, maxIter, r);
            int idx = zIter.Count - 1;
            if (maxIdx < idx)
                maxIdx = idx;
            xyIdx[i].Add(j, idx);
    for (int i = 0; i < w; i++)
        for (int j = 0; j < h; j++)
            int idx = xyIdx[i][j];
            double x = xMin + i * xStep;
            double y = yMin + j * yStep;
            ComplexNumber z = new ComplexNumber(x, y);
            bmp.SetPixel(w - i - 1, j, ComplexHeatMap(idx, 0, maxIdx, z, r));
    return bmp;
private static IList<ComplexNumber> SqPolyIteration(ComplexNumber z0, ComplexNumber c, int n, double r = 0)
    IList<ComplexNumber> res = new List<ComplexNumber>();
    for (int i = 0; i < n; i++)
        if (r > 0)
            if (res.Last().Mod > r)
        res.Add(res.Last() * res.Last() + c);
    return res;
private static double CalculateR(ComplexNumber c)
    return (1 + Math.Sqrt(1 + 4*c.Mod))/2;
public static Color ComplexHeatMap(decimal value, decimal min, decimal max, ComplexNumber z, double r)
    decimal val = (value - min) / (max - min);
    return Color.FromArgb(
        Convert.ToByte(255 * val),
        Convert.ToByte(255 * (1 - val)),
        Convert.ToByte(255 * (z.Mod / r > 1 ? 1 : z.Mod / r))

The Results

Here are a few 5000x5000 pixels pictures made with this function.

Now, let’s zoom the Julia set for the f(z) = z2 — 0.74543 + 0.11301i polynomial, the whole set is contained in the sphere of radius 1.50197192317588.

Julia set

Julia set

Julia set

Julia set

Julia set

Julia set

It seems that the accuracy limit is reached in this picture, and all red elements belong to the Julia set. But it’s far from it. Increase the maxIter parameter and get even better approximation quality. Anyway, you can go on like this to infinity. No kidding :-)

Julia set

Subscribe to Kukuruku Hub

Or subscribe with RSS


Thank you for the article! Love it!
Awesome fractals! Do you have the source code on github/bitbucket/etc.?

Read Next