Julia Set

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:

in which ○ is a function composition.

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

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:

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:

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).

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

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).

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:

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:

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

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:

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) = z² — 1 + 0.213i polynomial and build a trajectory for 100 iterations at

• z₀ = 0.3 + 0.2i (red)
  
• z₀ = 0.31 + 0.2i (blue)
  

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) = z² — 1 + 0.2i polynomial (we changed the c constant a bit) and build a trajectory for 100 iterations, when 

• z₀ = 0.3 + 0.2i (red)
  
• z₀ = 0.31 + 0.2i (blue)
  

Fatou, Julia and Mandelbrot Sets

• The Fatou set of the f(z) = z² + c polynomial is a subset of a set of complex numbers, for each point of which the behavior of a function under iteration is regular. That is to say, trajectory of points generated by iterations, will not change much when changing the initial conditions in the area of starting point. Denoted as F(f).
  
• The Julia set of the f(z) = z² + c polynomial is a subset of a set of complex numbers, for each point of which the behavior of a function under iteration is chaotic. That is to say, small changes in the area of starting point can much affect the trajectory. Denoted as J(f).
  
• The Mandelbrot set is a set of c parameters of the f(z) = z² + c polynomial, for which the Julia set is connected.
  

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.

• An Attractor is a set of numerical values toward which a system tends to evolve. In our case, it’s a subset of fixed points, to which sequences of the function iterations converge, independently of the initial conditions.
  
• An attractor's basin of attraction of point z of function f is a subset of points in the area of z (denoted as A(z)), in which any trajectory started in one of the points converges to point z.
  
• An attractor’s basin of attraction to infinity A(∞) is a set of z points, for which trajectories go to infinity. In other words, they do not converge to some fixed point and don’t oscillate between two or more periodic points (a fixed point is a periodic point with a period equal to one).
  

Let’s consider A(∞) in the context of our discussions for some f(z) = z² + 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:

The last theorem we need sounds the following way:

• Let there be given a quadratic polynomial of the f(z) = z² + c form
  
• Let’s denote 

• Then, the following condition will be satisfied for  point and n > 0:

• That is to say, if the module of the n-th iteration is greater that R, f function iterations diverge. Which is equivalent to the fact that point z₀ enters the attractor’s basin of attraction to infinity. This implies that z₀ is not a point of the filled Julia set.
  

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) = z²+ 1, threshold R is equal to the golden ratio о_0.

The Algorithm

1. Select с to define the f(z) = z² + с polynomial
   
2. Compute R for the defined f(z) = z² + с 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

/// 

/// Generate bmp with Julia set
/// 
/// parameter of square dynamics
/// width of bmp
/// height of bmp
/// max iterations of function
/// window in complex plane
/// window in complex plane
/// window in complex plane
/// window in complex plane
/// 
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> xyIdx = new Dictionary>();
    int maxIdx = 0;
    for (int i = 0; i ());
        for (int j = 0; j  zIter = SqPolyIteration(z, c, maxIter, r);
            int idx = zIter.Count - 1;
            if (maxIdx  SqPolyIteration(ComplexNumber z0, ComplexNumber c, int n, double r = 0)
{
    IList res = new List();
    res.Add(z0);
    for (int i = 0; i  0)
        {
            if (res.Last().Mod > r)
            {
                break;
            }
        }
        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(
        255,
        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.

• f(z) = z² — 0.8 + 0.156i
  
• f(z) = z² + 0.285 + 0.01i
  
• f(z) = z² — 0.0085 + 0.71i
  

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

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 :-)