Chapter 1: History of Fractal Mathematics¶

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
from numba import jit, prange
import warnings
warnings.filterwarnings('ignore')

# Set default figure style
plt.rcParams['figure.figsize'] = [10, 8]
plt.rcParams['figure.dpi'] = 100

def cantor_set(ax, x, y, length, depth, max_depth):
    """Draw the Cantor set recursively, showing all stages."""
    # Draw current segment
    ax.plot([x, x + length], [y, y], color='black', lw=3)

    if depth < max_depth:
        new_length = length / 3
        new_y = y - 0.15
        # Recurse on left and right thirds
        cantor_set(ax, x, new_y, new_length, depth + 1, max_depth)
        cantor_set(ax, x + 2 * new_length, new_y, new_length, depth + 1, max_depth)

# Visualize Cantor set construction
fig, ax = plt.subplots(figsize=(12, 6))
ax.set_title('Cantor Set Construction (6 iterations)', fontsize=14)
cantor_set(ax, 0, 1, 1, 0, 6)
ax.set_xlim(-0.05, 1.05)
ax.set_ylim(-0.1, 1.1)
ax.set_aspect('equal')
ax.axis('off')

# Add annotations
for i in range(7):
    ax.text(-0.08, 1 - i * 0.15, f'n={i}', fontsize=10, va='center')

# Calculate dimension
D_cantor = np.log(2) / np.log(3)
print(f"Cantor Set Fractal Dimension: D = log(2)/log(3) = {D_cantor:.4f}")
plt.tight_layout()
plt.show()

Cantor Set Fractal Dimension: D = log(2)/log(3) = 0.6309

def koch_curve(p1, p2, depth):
    """Generate Koch curve points between p1 and p2."""
    if depth == 0:
        return [p1, p2]

    # Calculate intermediate points
    dx = p2[0] - p1[0]
    dy = p2[1] - p1[1]

    # Points dividing the segment into thirds
    a = p1
    b = (p1[0] + dx/3, p1[1] + dy/3)
    d = (p1[0] + 2*dx/3, p1[1] + 2*dy/3)
    e = p2

    # Peak of the equilateral triangle
    # Rotate the middle third by 60 degrees
    cx = b[0] + dx/6 - dy * np.sqrt(3)/6
    cy = b[1] + dy/6 + dx * np.sqrt(3)/6
    c = (cx, cy)

    # Recurse on each segment
    points = []
    for seg in [(a, b), (b, c), (c, d), (d, e)]:
        pts = koch_curve(seg[0], seg[1], depth - 1)
        points.extend(pts[:-1])
    points.append(e)
    return points

def koch_snowflake(depth):
    """Generate a complete Koch snowflake."""
    # Start with equilateral triangle
    h = np.sqrt(3) / 2
    vertices = [(0, 0), (1, 0), (0.5, h)]

    all_points = []
    for i in range(3):
        p1 = vertices[i]
        p2 = vertices[(i + 1) % 3]
        points = koch_curve(p1, p2, depth)
        all_points.extend(points[:-1])
    all_points.append(all_points[0])  # Close the curve
    return all_points

# Visualize Koch curve evolution
fig, axes = plt.subplots(2, 3, figsize=(14, 8))
axes = axes.flatten()

for i, depth in enumerate(range(6)):
    points = koch_snowflake(depth)
    xs = [p[0] for p in points]
    ys = [p[1] for p in points]

    axes[i].plot(xs, ys, 'b-', lw=0.5 if depth > 3 else 1)
    axes[i].fill(xs, ys, alpha=0.3, color='lightblue')
    axes[i].set_title(f'Iteration {depth}', fontsize=12)
    axes[i].set_aspect('equal')
    axes[i].axis('off')

fig.suptitle('Koch Snowflake Construction', fontsize=14)

# Calculate dimension
D_koch = np.log(4) / np.log(3)
print(f"Koch Curve Fractal Dimension: D = log(4)/log(3) = {D_koch:.4f}")
print(f"Each iteration multiplies perimeter by 4/3 = {4/3:.4f}")
plt.tight_layout()
plt.show()

Koch Curve Fractal Dimension: D = log(4)/log(3) = 1.2619
Each iteration multiplies perimeter by 4/3 = 1.3333

def sierpinski_triangle(ax, vertices, depth, patches):
    """Generate Sierpinski triangle recursively."""
    if depth == 0:
        triangle = Polygon(vertices, closed=True)
        patches.append(triangle)
    else:
        # Calculate midpoints
        midpoints = [(vertices[i] + vertices[(i + 1) % 3]) / 2 for i in range(3)]
        # Recurse on three corner triangles
        sierpinski_triangle(ax, np.array([vertices[0], midpoints[0], midpoints[2]]), depth-1, patches)
        sierpinski_triangle(ax, np.array([vertices[1], midpoints[1], midpoints[0]]), depth-1, patches)
        sierpinski_triangle(ax, np.array([vertices[2], midpoints[2], midpoints[1]]), depth-1, patches)

# Visualize Sierpinski triangle evolution
fig, axes = plt.subplots(2, 3, figsize=(14, 10))
axes = axes.flatten()

vertices = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3)/2]])

for i, depth in enumerate(range(6)):
    patches = []
    sierpinski_triangle(axes[i], vertices, depth, patches)
    collection = PatchCollection(patches, facecolor='darkblue', edgecolor='black', linewidth=0.5)
    axes[i].add_collection(collection)
    axes[i].set_title(f'Iteration {depth} ({3**depth} triangles)', fontsize=11)
    axes[i].set_xlim(-0.1, 1.1)
    axes[i].set_ylim(-0.1, 1)
    axes[i].set_aspect('equal')
    axes[i].axis('off')

fig.suptitle('Sierpinski Triangle Construction', fontsize=14)

D_sierpinski = np.log(3) / np.log(2)
print(f"Sierpinski Triangle Dimension: D = log(3)/log(2) = {D_sierpinski:.4f}")
plt.tight_layout()
plt.show()

Sierpinski Triangle Dimension: D = log(3)/log(2) = 1.5850

def sierpinski_carpet(depth, size=243):
    """Generate Sierpinski carpet as a binary image."""
    # Use size that's a power of 3 for clean division
    carpet = np.ones((size, size), dtype=bool)

    def remove_center(x, y, s):
        """Recursively remove center squares."""
        if s < 3:
            return
        third = s // 3
        # Remove center square
        carpet[y + third:y + 2*third, x + third:x + 2*third] = False
        # Recurse on 8 surrounding squares
        for dy in range(3):
            for dx in range(3):
                if dx == 1 and dy == 1:
                    continue  # Skip center
                remove_center(x + dx*third, y + dy*third, third)

    for d in range(depth):
        step = size // (3 ** (d + 1))
        for y in range(0, size, step * 3):
            for x in range(0, size, step * 3):
                carpet[y + step:y + 2*step, x + step:x + 2*step] = False

    return carpet

# Visualize Sierpinski carpet evolution
fig, axes = plt.subplots(2, 3, figsize=(14, 10))
axes = axes.flatten()

for i, depth in enumerate(range(1, 7)):
    carpet = sierpinski_carpet(depth)
    axes[i].imshow(carpet, cmap='binary', interpolation='nearest')
    axes[i].set_title(f'Iteration {depth}', fontsize=11)
    axes[i].axis('off')

fig.suptitle('Sierpinski Carpet Construction', fontsize=14)

D_carpet = np.log(8) / np.log(3)
print(f"Sierpinski Carpet Dimension: D = log(8)/log(3) = {D_carpet:.4f}")
plt.tight_layout()
plt.show()

def weierstrass(x, a=0.5, b=7, n_terms=50):
    """Compute the Weierstrass function.

    Parameters:
    - a: amplitude decay (0 < a < 1)
    - b: frequency multiplier (odd integer)
    - n_terms: number of terms in the series
    """
    result = np.zeros_like(x)
    for n in range(n_terms):
        result += (a ** n) * np.cos((b ** n) * np.pi * x)
    return result

# Visualize Weierstrass function at different scales
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Different zoom levels
ranges = [(0, 2), (0, 0.5), (0, 0.1), (0, 0.02)]
titles = ['Full view [0, 2]', 'Zoom [0, 0.5]', 'Zoom [0, 0.1]', 'Zoom [0, 0.02]']

for ax, (x_min, x_max), title in zip(axes.flatten(), ranges, titles):
    x = np.linspace(x_min, x_max, 5000)
    y = weierstrass(x, a=0.5, b=7)
    ax.plot(x, y, 'b-', lw=0.5)
    ax.set_title(f'Weierstrass Function: {title}', fontsize=11)
    ax.set_xlabel('x')
    ax.set_ylabel('W(x)')
    ax.grid(True, alpha=0.3)

fig.suptitle('Weierstrass Function: Continuous but Nowhere Differentiable\n(a=0.5, b=7)', fontsize=14)
plt.tight_layout()
plt.show()

# The fractal dimension relates to parameters
a, b = 0.5, 7
D_weierstrass = 2 + np.log(a) / np.log(b)
print(f"Weierstrass function graph dimension: D ≈ 2 + log(a)/log(b) = {D_weierstrass:.4f}")

@jit(nopython=True, parallel=True)
def compute_julia(c_real, c_imag, x_min, x_max, y_min, y_max, width, height, max_iter):
    """Compute Julia set using escape-time algorithm."""
    result = np.zeros((height, width))

    for j in prange(height):
        for i in range(width):
            # Map pixel to complex plane
            x = x_min + (x_max - x_min) * i / width
            y = y_min + (y_max - y_min) * j / height

            # Iterate z -> z^2 + c
            zx, zy = x, y
            iteration = 0

            while zx*zx + zy*zy < 4 and iteration < max_iter:
                xtemp = zx*zx - zy*zy + c_real
                zy = 2*zx*zy + c_imag
                zx = xtemp
                iteration += 1

            result[j, i] = iteration

    return result

# Gallery of Julia sets with different c values
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.flatten()

# Famous Julia set parameters
c_values = [
    (-0.7, 0.27015, "Dendrite"),
    (-0.8, 0.156, "Rabbit"),
    (-0.4, 0.6, "Spiral"),
    (0.285, 0.01, "Siegel disk"),
    (-0.835, -0.2321, "Parabolic"),
    (-0.7269, 0.1889, "Dragon")
]

for ax, (c_re, c_im, name) in zip(axes, c_values):
    julia = compute_julia(c_re, c_im, -1.5, 1.5, -1.5, 1.5, 600, 600, 256)
    ax.imshow(julia, cmap='hot', extent=[-1.5, 1.5, -1.5, 1.5])
    ax.set_title(f'{name}\nc = {c_re} + {c_im}i', fontsize=10)
    ax.axis('off')

fig.suptitle('Julia Sets: Gallery of Different c Parameters', fontsize=14)
plt.tight_layout()
plt.show()

@jit(nopython=True, parallel=True)
def compute_mandelbrot(x_min, x_max, y_min, y_max, width, height, max_iter):
    """Compute Mandelbrot set using escape-time algorithm."""
    result = np.zeros((height, width))

    for j in prange(height):
        for i in range(width):
            # Map pixel to complex plane (c value)
            c_real = x_min + (x_max - x_min) * i / width
            c_imag = y_min + (y_max - y_min) * j / height

            # Iterate z -> z^2 + c starting from z = 0
            zx, zy = 0.0, 0.0
            iteration = 0

            while zx*zx + zy*zy < 4 and iteration < max_iter:
                xtemp = zx*zx - zy*zy + c_real
                zy = 2*zx*zy + c_imag
                zx = xtemp
                iteration += 1

            result[j, i] = iteration

    return result

# Full Mandelbrot set
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# Different zoom levels showing self-similarity
regions = [
    ((-2.5, 1.0, -1.25, 1.25), "Full Mandelbrot Set"),
    ((-0.75, -0.65, 0.1, 0.2), "Seahorse Valley"),
    ((-0.745, -0.735, 0.12, 0.13), "Deeper Zoom"),
    ((-0.748, -0.744, 0.098, 0.102), "Mini Mandelbrot")
]

for ax, ((x_min, x_max, y_min, y_max), title) in zip(axes.flatten(), regions):
    mandelbrot = compute_mandelbrot(x_min, x_max, y_min, y_max, 800, 800, 500)
    ax.imshow(mandelbrot, cmap='twilight_shifted', extent=[x_min, x_max, y_min, y_max])
    ax.set_title(title, fontsize=11)
    ax.set_xlabel('Re(c)')
    ax.set_ylabel('Im(c)')

fig.suptitle('Mandelbrot Set: Self-Similarity at Different Scales', fontsize=14)
plt.tight_layout()
plt.show()

print("The Mandelbrot set boundary has Hausdorff dimension = 2")
print("Small copies of the entire set appear at every scale.")

The Mandelbrot set boundary has Hausdorff dimension = 2
Small copies of the entire set appear at every scale.

# Calculate and display all dimensions
fractals = [
    ("Cantor Set", 2, 3, 1883),
    ("Koch Curve", 4, 3, 1904),
    ("Sierpinski Triangle", 3, 2, 1915),
    ("Sierpinski Carpet", 8, 3, 1916),
]

print("Self-Similar Fractal Dimensions: D = log(N) / log(1/r)")
print("=" * 60)
print(f"{'Fractal':<25} {'Year':>6} {'N':>4} {'1/r':>5} {'Dimension':>12}")
print("-" * 60)

for name, N, inv_r, year in fractals:
    D = np.log(N) / np.log(inv_r)
    print(f"{name:<25} {year:>6} {N:>4} {inv_r:>5} {D:>12.4f}")

print("=" * 60)

Self-Similar Fractal Dimensions: D = log(N) / log(1/r)
============================================================
Fractal                     Year    N   1/r    Dimension
------------------------------------------------------------
Cantor Set                  1883    2     3       0.6309
Koch Curve                  1904    4     3       1.2619
Sierpinski Triangle         1915    3     2       1.5850
Sierpinski Carpet           1916    8     3       1.8928
============================================================