The Koch snowflake

How can a shape that fits inside a circle have an endless boundary? The Koch snowflake starts from an equilateral triangle. In every step, the middle third of each edge is replaced by two sides of a small equilateral triangle pointing outward. Every edge becomes four edges a third as long, so the perimeter grows by \(4/3\) per step, without end. The area grows too, but by less each time, and tends to \(8/5\) of the triangle’s area \(A_0\). After \(n\) steps,

\[A_n = A_0 \left(1 + \tfrac{3}{5} \left(1 - \left(\tfrac{4}{9}\right)^n\right)\right), \qquad P_n = 3 s \left(\tfrac{4}{3}\right)^n.\]

The example builds six steps as polylines and lets planex2d compute the area and the perimeter of each, which are checked against the formulas. It also checks that every step stays inside the circle through the triangle’s corners.

This page is part of a collection of geometry and design examples using planex2d.

The Koch snowflake after five steps, the first steps drawn inside it

[1]:
import math

import numpy as np
import planex as px

The steps

The triangle has the side \(s = 100\) and runs counterclockwise, so outward is to the right of each edge. A step turns the middle third of every edge by 60° to the right and adds the corner of the small triangle. NumPy does this for all edges at once:

[3]:
s = 100
triangle = np.array([(0, 0), (s, 0), (s / 2, s * math.sqrt(3) / 2)])


def koch(points):
    """every edge's middle third replaced by two sides of an equilateral triangle, pointing outward"""
    third = (np.roll(points, -1, axis=0) - points) / 3
    c, si = math.cos(math.radians(-60)), math.sin(math.radians(-60))
    turned = np.column_stack([c * third[:, 0] - si * third[:, 1], si * third[:, 0] + c * third[:, 1]])
    return np.stack([points, points + third, points + third + turned, points + 2 * third], axis=1).reshape(-1, 2)

Area and perimeter

Each step is a closed polyline, and planex2d’s face inside it gives the area and the length of its outline (Figure 1):

[4]:
A0 = math.sqrt(3) / 4 * s**2
points = triangle
flakes = {}
for n in range(7):
    flakes[n] = px.FromPolyline(points, closed=True).Face()
    print("n", n, "| corners", len(points), "| area / A0", flakes[n].mass / A0, "formula",
          1 + 3 / 5 * (1 - (4 / 9)**n), "| perimeter", flakes[n].wires[0].mass, "3 s (4/3)^n", 3 * s * (4 / 3)**n)
    points = koch(points)
# px.Plot(flakes[6])
n 0 | corners 3 | area / A0 1.0 formula 1.0 | perimeter 300.0 3 s (4/3)^n 300.0
n 1 | corners 12 | area / A0 1.333333333333333 formula 1.3333333333333333 | perimeter 399.9999999999999 3 s (4/3)^n 400.0
n 2 | corners 48 | area / A0 1.4814814814814814 formula 1.4814814814814814 | perimeter 533.3333333333331 3 s (4/3)^n 533.3333333333333
n 3 | corners 192 | area / A0 1.5473251028806572 formula 1.5473251028806585 | perimeter 711.1111111111097 3 s (4/3)^n 711.111111111111
n 4 | corners 768 | area / A0 1.5765889346136224 formula 1.576588934613626 | perimeter 948.1481481481372 3 s (4/3)^n 948.1481481481479
n 5 | corners 3072 | area / A0 1.5895950820504894 formula 1.5895950820505003 | perimeter 1264.1975308642545 3 s (4/3)^n 1264.197530864197
n 6 | corners 12288 | area / A0 1.5953755920224244 formula 1.5953755920224446 | perimeter 1685.5967078190215 3 s (4/3)^n 1685.5967078189294
After 1 step: 12 edges
(a) After 1 step: 12 edges
After 2 steps: 48 edges
(b) After 2 steps: 48 edges
After 3 steps: 192 edges
(c) After 3 steps: 192 edges
After 6 steps: 12288 edges
(d) After 6 steps: 12288 edges
Figure 1: The Koch snowflake step by step

Inside a circle

The corners of the small triangles never reach past the circle through the triangle’s corners, so the whole snowflake stays inside it, and the Boolean difference of each step and the disc is empty. Its area stays below \(8/5\) of the triangle’s while its perimeter grows without end, so the isoperimetric quotient \(4 \pi A / P^2\), which is 1 for the circle, falls toward zero:

[6]:
center = (s / 2, s * math.sqrt(3) / 6)
disc = px.Circle(center, s / math.sqrt(3)).Face()
for n, flake in flakes.items():
    print("n", n, "| outside the disc", (flake - disc).mass, "| 4πA/P²", 4 * math.pi * flake.mass / flake.wires[0].mass**2)
n 0 | outside the disc 0.0 | 4πA/P² 0.6045997880780726
n 1 | outside the disc 0.0 | 4πA/P² 0.4534498410585546
n 2 | outside the disc 0.0 | 4πA/P² 0.2834061506615967
n 3 | outside the disc 0.0 | 4πA/P² 0.1665011135136885
n 4 | outside the disc 0.0 | 4πA/P² 0.09542816479308643
n 5 | outside the disc 0.0 | 4πA/P² 0.05412116480651349
n 6 | outside the disc 0.0 | 4πA/P² 0.03055386073126528

See The isoperimetric inequality for the quotient and The coastline paradox for a real boundary whose length depends on the ruler.