In [1]:
from matplotlib import pyplot as plt
from matplotlib import collections as mc
import math
%matplotlib inline
A Koch Snowflake is a fractal that has been known for over 100 years (see the Wikipedia article for history).
The shaped is formed by starting from a triangle. For each line segment, remove the middle third and replace it by two equal pieces that form a triangle.
To program this in python, we simply need a function that turns a line segment into four shorter segments. We'll use a pair of tuples to represent a line segment: $((x_a,y_a), (x_b,y_b))$. The current shape will be a list of segments.
All we need is a function that takes one line segment and expands it into four smaller segments. We'll call the original segment ae, and create new points b, c, d, so that the four segemnts are ab, bc, cd, and de. First, let's work this out by hand, then we'll make the function.
In [2]:
a = (0.0, 0.0)
e = (1.0, 0.0)
ae = (a,e)
It's helpful to be able to make a plot of a list of segments.
In [3]:
def plot_segments(segments):
fig, ax = plt.subplots()
lines = mc.LineCollection(segments)
ax.add_collection(lines)
ax.margins(0.2)
ax.set_aspect('equal')
ax.autoscale()
return ax
In [4]:
plot_segments([ae]);
Next we figure out formulas for points b, c, and d. Points b and d are easy, because they are 1/3 and 2/3 of the way along the segment ae.
In [5]:
b = ((2*a[0]+e[0]/3, (2*a[1]+e[1])/3))
d = ((a[0]+2*e[0]/3, (a[1]+2*e[1])/3))
Point c is trickier, because it doesn't lie directly on the line segment. It is the vertex of an equilateral triangle with side length |ae|/3. To get to point c, find the midpoint of ae, then go out perpendicularly a distance $\sqrt{3}/6$. To move perpendicularly, we use that trick that the point (-y, x) is rotated 90° CCW from (x, y).
In [6]:
k = math.sqrt(3)/6
c = ((a[0]+e[0])/2 - k * (e[1]-a[1]), (a[1]+e[1])/2 + k *(e[0]-a[0]))
In [7]:
plt.gcf().clear()
plot_segments([(a,b), (b,c), (c,d), (d,e)]);
Now we make this into a function.
In [8]:
def f(seg):
a = seg[0]
e = seg[1]
b = ((2*a[0]+e[0])/3, (2*a[1]+e[1])/3)
d = ((a[0]+2*e[0])/3, (a[1]+2*e[1])/3)
k = math.sqrt(3)/6
c = ((a[0]+e[0])/2 - k * (e[1]-a[1]), (a[1]+e[1])/2 + k *(e[0]-a[0]))
return [(a,b), (b,c), (c,d), (d,e)]
We'll test this function on some different line segments.
In [9]:
plot_segments(f(((0,0),(1,0))));
In [10]:
plot_segments(f(((0,0),(0,1))));
In [11]:
plot_segments(f(((2, 3), (2 + math.cos(math.pi/3), 3 + math.sin(math.pi/3)))));
Finally, we make a function to apply f to every segment in a list. We use some elegant notation called a “list comprehension“ here.
In [12]:
def recurse(segments):
return [x for s in segments for x in f(s)]
In [13]:
recurse([(a,e)])
Out[13]:
In [14]:
plot_segments(recurse([(a,e)]));
In [15]:
segements = [(a,e)]
for i in range(2):
segements = recurse(segements)
plot_segments(segements);
Finally, we'll make the full snowflake by starting from an equilateral triangle.
In [16]:
def snowflake(n):
p = -math.cos(math.pi/6), math.sin(math.pi/6)
q = math.cos(math.pi/6), math.sin(5*math.pi/6)
r = 0.0, -1.0
segments = [(p,q), (q,r), (r,p)]
for i in range(n):
segments = recurse(segments)
plot_segments(segments)
In [17]:
snowflake(0)
In [18]:
snowflake(1)
In [19]:
snowflake(2)
In [20]:
snowflake(6)
In [21]:
[(n, 3*(4/3)**n) for n in range(11)]
Out[21]:
Note that the true fractal curve, with $n\rightarrow\infty$), has an infinite length!
For the total area of the fractal, consider the additional area at each iteration. For simplicity, we'll measure area relative to the starting triangle:
To continue this, we need to be more systematic. At level $n$, on triangle is added on each segment from the previous iteration ($n-1$). The number of segments at level $n$ is $3\cdot 4^n$, so the number of triangles added at level $n$ is $3\cdot 4^{n-1}$. Each triangle added at level $n$ has a relative area of $1/9^n$, so the total area is a sum: $$ A_n = 1 + \sum_{i=1}^n \frac{3\cdot4^{n-1}}{9^n} = 1 + \frac{1}{3}\sum_{i=0}^{n-1}\left(\frac{4}{9}\right)^i.$$
Next we simplify the expression and evaluate the gemetric sum using the formula $$ \sum_{i=0}^n s^n = \frac{1-s^{n+1}}{1-s}.$$ We find $$A_n = 1 + \frac{1}{3}\sum_{i=0}^{n-1} \left(\frac{4}{9}\right)^i = 1 + \frac{1}{3}\cdot\frac{1-(4/9)^n}{1-(4/9)} = \frac{8}{5} - \frac{3}{5}\left(\frac{4}{9}\right)^n$$
In [22]:
[(n, 8/5 - 3/5*(4/9)**n) for n in range(11)]
Out[22]:
The series exponentially converges to the limit $8/5 = 1.6$.