Mandelbrot fractal area approximation by Monte Carlo sampling

monte_carlo_mandelbrot_area.py

"""
Approximates area of Mandelbrot fractal by Monte Carlo sampling.

The definition of the Mandelbrot Set is the set of values c \in C where

  z_{n+1} = z_n^2 + c

remains bounded. Most members will be within (-2, 2) x (-2i, 2i) so
here this is used as an artificial bounds to simplify the approximation.

"""
import numpy as np
from PIL import Image, ImageDraw

# Number of iterations after which z_n is considered bounded.
iters_bounded = 200
N = 100000
width = 400
height = 400

im = Image.new("RGB", (width, height), color="black")
draw = ImageDraw.Draw(im)
indicator_count = 0

# Sample N complex coordinates to find the ratio that is within
# Mandelbrot fractal.
for px in range(N):
    c = np.complex(np.random.uniform(-2, 2), np.random.uniform(-2, 2))
    z = 0
    j = 0
    while j < iters_bounded and np.abs(z) < 2:
        z = z * z + c
        j += 1
    if j == iters_bounded:
        indicator_count += 1
        x = width / 2 + np.real(c) * width / 4
        y = height / 2 + np.imag(c) * height / 4
        draw.point((x, y), fill="white")

# Area of Mandelbrot is (-2, 2) x (-2, 2) = 16 and indicator_count / N is
# the fraction of samples (0, 1] that was within Mandelbrot fractal.
approx = 16 * indicator_count / N
print("Area is approximately", approx)

# Show samples visually.
im.show()