Complex Dynamics

Julia Set — c = −0.7 + 0.27i

The Julia set for a specific complex constant, computed by iterating z² + c until escape.

The Iteration

Pick any point (x, y) in the plane and treat it as a complex number:

z = x + yi

Now repeatedly apply this rule:

z_{n+1} = z_n² + c

where c = -0.7 + 0.27i is a fixed complex constant. The question is: does this sequence of values escape to infinity, or stay bounded forever?

Stage 0 shows the pixel grid before we’ve computed anything. Every pixel is a candidate — we haven’t yet determined which ones escape and which ones don’t.

Escape-Time Coloring

We iterate each pixel up to 256 times. If the magnitude |z| ever exceeds 2, the sequence will definitely escape to infinity. The escape time — how many iterations it took — becomes the color.

def iterate(x, y, c, max_iter=256):
    z = complex(x, y)
    for i in range(max_iter):
        if abs(z) > 2:
            return i      # escaped at iteration i
        z = z*z + c
    return max_iter       # bounded — part of the Julia set

Pixels that escape quickly get bright colors. Pixels that take longer are darker. Pixels that never escape (they’re in the Julia set itself, or its interior) are colored black.

The color bands you see are level curves — all the points that escaped at the same iteration count.

The Boundary

The Julia set itself is the boundary between the two regions: the set of all points z where the sequence is neither clearly escaping nor clearly bounded. It’s a fractal curve.

For c = -0.7 + 0.27i, the Julia set is a connected, intricate curve — this particular value of c produces a “dendrite” shape, dense with filaments.

The key insight: the shape of the Julia set is entirely determined by c. Change c by even a tiny amount and you get a completely different fractal. There is a whole universe of Julia sets — one for every point in the complex plane.

import numpy as np

def julia_grid(width, height, c, max_iter=256):
    x = np.linspace(-1.5, 1.5, width)
    y = np.linspace(-1.5, 1.5, height)
    Z = x[np.newaxis, :] + 1j * y[:, np.newaxis]

    count = np.zeros(Z.shape, dtype=int)
    mask  = np.ones(Z.shape, dtype=bool)   # True = still iterating

    for i in range(max_iter):
        Z[mask] = Z[mask]**2 + c
        escaped = mask & (np.abs(Z) > 2)
        count[escaped] = i
        mask[escaped] = False

    return count

The key to performance is operating on the entire grid simultaneously using NumPy arrays rather than looping pixel by pixel. This vectorised approach renders a 600×600 grid roughly 100× faster than a pure Python loop.