How to Make Fractals in Python

There is something vertiginous about writing twenty lines of Python and watching your terminal produce a shape with infinite boundary. That is what fractal programming offers: a direct line between elementary arithmetic and one of the most complex geometric objects in mathematics. This guide walks you through four distinct Python techniques — from a minimal Mandelbrot renderer to an animated Barnsley fern — explaining both the code and the mathematics that drives each one.

What libraries do you need to make fractals in Python?

Before writing a single loop, it helps to know which tool fits which fractal. The table below maps each fractal family to its ideal Python stack:

All four libraries ship with a standard Python 3 installation or are installable in seconds via pip install numpy matplotlib. No exotic dependencies are required.

Can you make the Mandelbrot set in Python?

Yes — and it takes fewer lines than most people expect. The Mandelbrot set is defined by the iteration z_n+1 = z_n² + c, where both z and c are complex numbers and z starts at zero. For each point c on the complex plane, you iterate the formula repeatedly: if |z| stays bounded (never exceeds 2), the point is in the set and colored black; if it escapes, you color it by how quickly it left. The result is the iconic shape first visualized on an IBM supercomputer in 1980.

The naive pure-Python implementation uses a pixel-by-pixel loop:

import numpy as np
import matplotlib.pyplot as plt

width, height, max_iter = 800, 600, 256
x = np.linspace(-2.5, 1.0, width)
y = np.linspace(-1.25, 1.25, height)
C = x[np.newaxis, :] + 1j * y[:, np.newaxis]  # complex grid

Z = np.zeros_like(C)
M = np.zeros(C.shape, dtype=int)

for i in range(max_iter):
    mask = np.abs(Z) <= 2
    Z[mask] = Z[mask] ** 2 + C[mask]
    M[mask] += 1

plt.figure(figsize=(10, 7))
plt.imshow(M, extent=(-2.5, 1.0, -1.25, 1.25), cmap="inferno", origin="lower")
plt.title("Mandelbrot Set")
plt.colorbar(label="Escape iteration")
plt.tight_layout()
plt.show()

The key insight in this vectorized version — drawn from the authoritative Real Python tutorial and the official NumPy fractal tutorial — is that the mask array prevents us from updating points that have already escaped: once a pixel's z exceeds 2 it is permanently excluded from future iterations. This single boolean mask turns a loop that would otherwise touch 480,000 complex-number multiplications per iteration into a shrinking operation that speeds up automatically as more points escape. At max_iter=256 on modern hardware the full image renders in under two seconds.

Color choice matters more than it might seem. The "inferno" colormap assigns dark purples to points that escape quickly (far outside the set) and bright yellows to points that linger longest at the boundary — visually emphasizing the intricate filaments that define the Mandelbrot set's fractal boundary, whose fractal dimension is exactly 2.

How do you draw a fractal tree and Koch snowflake in Python?

Geometric fractals like the Koch snowflake and fractal trees are built on the mathematical concept of self-similarity: each piece of the structure is a scaled copy of the whole. Python's built-in turtle module — which ships with every standard Python 3 installation — mirrors this structure naturally, because turtle graphics are themselves recursive: to draw a branch, you draw a shorter branch at an angle, then another, then another, until the branches become too small to see.

Here is a classic recursive fractal tree:

import turtle

def draw_branch(t, length, angle=25, scale=0.7):
    if length < 5:
        return
    t.forward(length)
    t.right(angle)
    draw_branch(t, length * scale, angle, scale)
    t.left(angle * 2)
    draw_branch(t, length * scale, angle, scale)
    t.right(angle)
    t.backward(length)

screen = turtle.Screen()
screen.bgcolor("black")
t = turtle.Turtle()
t.color("lime")
t.speed(0)
t.left(90)
t.up()
t.goto(0, -200)
t.down()
draw_branch(t, 120)
turtle.done()

The critical structural element is the pair of recursive calls separated by heading adjustments: t.right(angle) before the first recursive call turns the turtle to draw the right sub-branch; t.left(angle * 2) overshoots back past center to draw the left sub-branch; then t.right(angle) restores the original heading. The t.backward(length) at the base of the function is equally important — it walks the turtle back to the branch's starting point so higher levels of the recursion can continue correctly. Without this backtracking, successive branches would drift off into open space.

For the Koch snowflake, the same recursive logic applies, but the "branching" rule is different: each straight segment is replaced by four segments (the middle third of the original is replaced by two sides of an equilateral triangle):

def koch_curve(t, length, depth):
    if depth == 0:
        t.forward(length)
        return
    length /= 3
    koch_curve(t, length, depth - 1)
    t.left(60)
    koch_curve(t, length, depth - 1)
    t.right(120)
    koch_curve(t, length, depth - 1)
    t.left(60)
    koch_curve(t, length, depth - 1)

A depth of 5 produces a visually rich snowflake; depth 7 or higher starts to reveal the paradox the Koch snowflake illustrates: each iteration multiplies the perimeter by 4/3, so after n iterations the perimeter is (4/3)ⁿ times the original — growing without bound even as the enclosed area converges to a finite value. This geometric fact connects directly to the fractal dimension of the Koch curve, which is log(4)/log(3) ≈ 1.262 — strictly between the topological dimension of a line (1) and a filled plane (2).

The official Matplotlib blog demonstrates how to extend this approach to animated fractals, successively revealing each iteration of the Koch snowflake or Mandelbrot set using FuncAnimation.

How do you generate a Barnsley fern in Python?

The Barnsley fern, described by mathematician Michael Barnsley in his 1988 textbook Fractals Everywhere, is a fractal that genuinely looks like a plant — and it is generated by a completely different mechanism than either escape-time or recursive geometric fractals. It uses an Iterated Function System (IFS): four affine transformations are applied randomly with fixed probabilities, each one "pushing" a point into a different region of the final fern shape. Plot enough points (typically 100,000 or more) and the emergent structure is unmistakable.

import random
import matplotlib.pyplot as plt

def barnsley_fern(n=200_000):
    x, y = 0.0, 0.0
    xs, ys = [], []
    for _ in range(n):
        r = random.random()
        if r < 0.01:          # stem
            x, y = 0.0, 0.16 * y
        elif r < 0.86:        # main leaflet
            x, y = 0.85 * x + 0.04 * y, -0.04 * x + 0.85 * y + 1.6
        elif r < 0.93:        # right leaflet
            x, y = 0.20 * x - 0.26 * y,  0.23 * x + 0.22 * y + 1.6
        else:                  # left leaflet
            x, y = -0.15 * x + 0.28 * y, 0.26 * x + 0.24 * y + 0.44
        xs.append(x)
        ys.append(y)
    return xs, ys

xs, ys = barnsley_fern()
plt.figure(figsize=(6, 10), facecolor="black")
plt.scatter(xs, ys, s=0.1, c=ys, cmap="Greens", linewidths=0)
plt.axis("off")
plt.title("Barnsley Fern", color="white")
plt.tight_layout()
plt.show()

The four transformation matrices encode the biological structure of a fern: the first (probability 1%) maps everything to the stem; the second (85%) handles the main compound leaf axis; the third and fourth (7% and 7%) generate the sub-leaflets. The magic of the IFS is that the attractor — the shape the random walk converges toward — is entirely determined by the matrices, not by the starting point. You can begin at any (x, y) and after a few hundred "warm-up" iterations the trajectory locks onto the fern. This is a concrete demonstration of the self-similarity principle: each leaflet is a miniaturized, transformed copy of the whole fern.

Coloring the scatter plot by the y-coordinate (height) using a green colormap produces depth cues that reinforce the biological illusion. At 200,000 points the image is crisp; at 500,000 points individual leaflet veins become visible.

How do you speed up Python fractal rendering with NumPy?

Pure Python loops over pixel grids are slow — a 1,000 × 1,000 Mandelbrot image iterated 256 times requires roughly 256 million arithmetic operations. NumPy's vectorized operations run in compiled C code under the hood, typically achieving 50–200× speedups over pure Python on the same calculation.

The core optimization is replacing per-pixel loops with array-wide operations. Instead of:

# Slow: pure Python loop
for i in range(height):
    for j in range(width):
        result[i, j] = mandelbrot(C[i, j], max_iter)

you operate on the entire 2D complex array at once, using a boolean mask to exclude already-escaped points:

# Fast: NumPy vectorized
Z = np.zeros_like(C)
M = np.zeros(C.shape, dtype=int)
for i in range(max_iter):
    mask = np.abs(Z) <= 2          # points still inside
    Z[mask] = Z[mask]**2 + C[mask]  # iterate only bounded points
    M[mask] += 1                    # count iterations

For even greater performance, the SciPy lecture notes demonstrate how to use Numba's @jit decorator to JIT-compile the inner loop to native machine code, often achieving speeds within 5× of hand-written C. For very high-resolution renders (4K and above), the multiprocessing module can distribute rows of the image across CPU cores, since each row is computed independently.

A practical reference point: on a modern laptop, a 2,000 × 2,000 Mandelbrot at 512 iterations takes roughly 40 seconds in pure Python, 1.5 seconds with vectorized NumPy, and under 0.3 seconds with Numba JIT. For exploratory work, NumPy is almost always sufficient.

What are the best Python packages for fractal visualization beyond Matplotlib?

Matplotlib is the workhorse for static images, but the Python fractal ecosystem extends considerably beyond it:

• HoloViews + Datashader — handles billions of IFS points without memory exhaustion, ideal for ultra-high-density Barnsley fern or flame fractal renders.
• Pillow (PIL) — produces direct pixel-level image output without a display window; useful for saving fractal images to disk programmatically. The Real Python Mandelbrot tutorial demonstrates Pillow rendering alongside Matplotlib.
• Pygame — enables real-time interactive fractal explorers with keyboard zoom controls; the rendering loop maps naturally onto the escape-time algorithm.
• Manim — the mathematical animation library (behind the 3Blue1Brown YouTube channel) can animate fractal construction step by step, ideal for educational visualizations.
• mpmath — Python's arbitrary-precision mathematics library, essential when zooming deep into the Mandelbrot set where standard 64-bit floats lose resolution (typically below magnifications of 10¹⁵).

For three-dimensional fractals — Mandelbulbs, Menger sponges, Apollonian gaskets — the dedicated software Mandelbulb 3D and Mandelbulber are purpose-built and faster than anything achievable in pure Python; but Python scripts can still drive them programmatically via their command-line interfaces.

If you want to share your fractals as interactive web visualizations, Plotly's imshow function renders escape-time fractals as zoomable, pannable HTML canvases directly exportable to a web page — no server required.