Koch Fractals

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

Koch curves are particularly simple fractals which enclose a finite area within a continuous, nonintersecting curve of infinite length.  A Koch curve is generated by the following self-similarity transformation.  Begin with a regular polygon composed of sides of equal length (Subscript[ℓ, 0]).  Centered on each segment of length (Subscript[ℓ, 0]) build an equilateral triangle with sides of length (Subscript[ℓ, 1]) = (Subscript[ℓ, 0])/3 using the middle third of the original segment as a base.  Thus, after one iteration a simple line segment becomes

Then build triangles with sides of length (Subscript[ℓ, 2]) = (Subscript[ℓ, 1])/3 upon the middle third of each segment of length (Subscript[ℓ, 1]), such that the preceding figure becomes

This transformation is repeated indefinitely and the Koch curve is the border of the resulting figure.

For example, the Koch snowflake based upon an equilateral triangle is shown below, where here the number of iterations is n = 4.  The figure on the left shows the Koch curve, which is only the outer border, while the figure on the right shows the intermediate steps to illustrate the construction algorithm.  In addition, we also show the circle that circumscribes the original triangle to illustrate the finiteness of the area.  However, one can show that in the limit n -> ∞ the length of the Koch curve is infinite despite the fact that it encloses a finite area.

Implementation

Initialization

ClearAll["Global`*"];
Off[General::spell,  General::spell1]

Needs["Utilities`Notation`"]

Symbolize[(Subscript[n, 0])]; Symbolize[(Subscript[ℓ, 0])]; Symbolize[(Subscript[L, 0])];Symbolize[(Subscript[A, 0])]

Basic algorithm

The vertices of the triangle built on a line segment

To construct an equilateral triangle upon a line segment, we must calculate the base points and the vertex point.  The simplest way to do this is with vectors.  Let (Overscript[r, ⇀]) = {x, y} represent a two-dimensional vector and let  (Subscript[(Overscript[r, ⇀]), 1]) = {x, y} and  (Subscript[(Overscript[r, ⇀]), 2]) = {x, y} represent the end points of a line segment.  The base points of the triangle are then \(TraditionalForm
      \`r\&⇀\_a = 
        r\&⇀\_1 + 
          \(1\/3\) \((r\&⇀\_2 - r\&⇀\_1)\)\) and  (Subscript[(Overscript[r, ⇀]), b]) = 
              (Subscript[(Overscript[r, ⇀]), 1]) + 
                (2/3) ((Subscript[(Overscript[r, ⇀]), 2]) - (Subscript[(Overscript[r, ⇀]), 1])).  To compute the vertex, we must move along a perpendicular to the line segment a distance equal to (1/(Sqrt[12])) its length.  Recognizing that the vector {y, -x} is perpendicular to {x, y}, we can compute the vertex point as

\(TraditionalForm\`\(r\& ⇀ \)\_c = \(1\/2\) \((\(r\& ⇀ \)\_1 + \(r\& ⇀ \)\_2)\) + {y\_2 - y\_1, x\_1 - x\_2}\/\@12\)

These points are evaluated using the following functions.

Clear[perp];
perp[
          vector : {x_, y_}] := {y, -x}/(Sqrt[(x^2) + (y^2)]);
perp[pt1 : {x1_, y1_}, pt2 : {x2_, y2_}] := perp[pt1 - pt2]

newpt1[pt1 : {x1_, y1_}, pt2 : {x2_, y2_}] := 
        pt1 + (1/3) (pt2 - pt1);

newpt2[pt1 : {x1_, y1_}, 
          pt2 : {x2_, y2_}] := (1/2) (pt1 + pt2) + (1/(Sqrt[12])) 
            perp[pt1, 
              pt2] (Sqrt[(pt2 - pt1) . (pt2 - pt1)]);

newpt2[{x1, y1}, {x2, y2}] // Simplify

newpt3[pt1 : {x1_, y1_}, pt2 : {x2_, y2_}] := 
        pt1 + (2/3) (pt2 - pt1);

Replacement rule for line segments

The following rule replaces a line segment by the four segments required to place a triangle atop the middle third of the original segment.

Clear[KochRule];
KochRule = 
        Line[{pt1 : {x1_, y1_}, 
              pt2 : {x2_, y2_}}] :> {Line[{pt1, 
                newpt1[pt1, pt2]}], 
            Line[{newpt1[pt1, pt2], newpt2[pt1, pt2]}], 
            Line[{newpt2[pt1, pt2], newpt3[pt1, pt2]}], 
            Line[{newpt3[pt1, pt2], pt2}]};

Clear[KochSplit];
KochSplit[x_] := x /. KochRule

Basic function which displays Koch transformation of polygon

Clear[BasicDisplayKoch];
BasicDisplayKoch[polygon : {{_, _} .. }, n_, opts___Rule] := 
      Show[Graphics[Nest[KochSplit, Line /@ Partition[polygon, 2, 1], n]], 
        opts, AspectRatio -> Automatic]

The basic algorithm is illustrated by the following iteration sequence applied to a simple line segment.

Table[
        BasicDisplayKoch[{{0. , 0. }, {1. , 0. }}, n], {n, 1, 
          4}];

More elaborate display function

Now that we have an algorithm which produces and displays Koch curves, we can write a more elaborate function which permits several variations to be chosen.  We use Options to establish default options for our function and the FilterOptions package to distinguish between our options and those applicable to Graphics.  We also provide usage messages for the function and its options.

Needs["Utilities`FilterOptions`"];

Clear[DisplayKoch];
\(\(DisplayKoch::usage = "DisplayKoch[polygon,n,opts] produces a Koch curve of order n based upon the list of points in polygon. By default, the triangle rule is applied in the outward direction assuming that the polygon is defined counterclockwise.

  Options: 

    KochType->"sieve" applies the triangle-building rule in reverse.

    ShowIterations->True can be used to display intermediate iterations.

    CloseFigure->True connects the last point to the first to close the figure.";\)
\)
\(\(ShowIterations::usage = "a True|False option for DisplayKoch.";\)
\)
\(\(KochType::usage = "an option for DisplayKoch: 

  "snowflake" produces outward triangles 

  "sieve" produces inward triangles";\)
\)
\(\(CloseFigure::usage = "a True|False option for DisplayKoch that can be used when the list of points is not a closed polygon.";\)
\)
\(\(Options[DisplayKoch] = {ShowIterations -> False, 
          KochType -> "snowflake", 
          CloseFigure -> False};\)
\)
\(DisplayKoch[polygon : {{_, _} .. }, n_, opts___Rule] := 
      Module[{curve, type, showIterations, close, NestFunction, optKoch, 
          optPlot}, 
	
 (*\ 
          establish\ options\ *) 

{type, showIterations, 
            close} = \({KochType, ShowIterations, CloseFigure} /. {opts}\) /. 
            Options[DisplayKoch]; 
		
 (*\ 
          establish\ initial\ line\ *) 


        curve = \ If[type == "sieve", Reverse[polygon], polygon]; 

        If[close, AppendTo[curve, First[curve]]]; 

        curve = Line /@ Partition[curve, 2, 1]; 
	
 (*\ 
          apply\ Koch\ algorithm\ and\ produce\ figure\ *) 

NestFunction = If[showIterations, NestList, Nest]; 

        optPlot = FilterOptions[Graphics, opts]; 

        Show[Graphics[NestFunction[KochSplit, curve, n]], optPlot, 
          AspectRatio -> Automatic]]\)

DisplayKoch[{{0. , 0. }, {0.5, ((Sqrt[3])/2.) }, {1. , 0. }}, 4, 
        ShowIterations -> False, KochType -> "sieve", 
        CloseFigure -> True];

Typical Koch curves

Simple figures based upon equilateral and square

The following specify an equilateral triangle and a square in a counterclockwise direction of traversal.

triangle = {{0. , 0. }, {0.5, ((Sqrt[3])/2.) }, {1. , 0. }, {0. , 
            0. }};
square = {{0. , 0. }, {1. , 0. }, {1. , 1. }, {0. , 1. }, {0. , 
            0. }};

The effect of outward versus inward triangulation are compared below for a (4^th) order Koch curve based upon the equilateral triangle.  The first pair of figures shows the intermediate iterations while the second pair shows the final result.

Show[
        GraphicsArray[{DisplayKoch[triangle, 4, ShowIterations -> True, 
              DisplayFunction -> Identity], 
            DisplayKoch[triangle, 4, ShowIterations -> True, 
              KochType -> "sieve", 
              DisplayFunction -> Identity]}], 
        DisplayFunction -> $DisplayFunction];

Show[
        GraphicsArray[{DisplayKoch[triangle, 4, 
              DisplayFunction -> Identity], 
            DisplayKoch[triangle, 4, KochType -> "sieve", 
              DisplayFunction -> Identity]}], 
        DisplayFunction -> $DisplayFunction];

Similar figures based upon a square are shown below.

Show[
        GraphicsArray[{DisplayKoch[square, 4, 
              DisplayFunction -> Identity], 
            DisplayKoch[square, 4, KochType -> "sieve", 
              DisplayFunction -> Identity]}], 
        DisplayFunction -> $DisplayFunction];

Figures based upon higher-order polygons

Needs["Geometry`Polytopes`"]

DisplayKoch[Vertices[Dodecagon], 4, ShowIterations -> False, 
        KochType -> "snowflake", 
        CloseFigure -> True];

Fractal properties

Infinite length encloses finite area

Consider a regular polygon with (Subscript[n, 0]) sides of length (Subscript[ℓ, 0]) enclosing area (Subscript[A, 0]).  Each iteration multiplies the number of line segments by 4 and reduces the length of each segment by a factor of (1/3).  Thus, the length after n iterations becomes (Subscript[L, n]) = ((((4/3)))^n) (Subscript[L, 0]) where (Subscript[L, 0]) = (Subscript[n, 0]) (Subscript[ℓ, 0]) and clearly diverges as n -> ∞.  Nevertheless, the following comparison between the circle that circumscribes the original triangle and a Koch triangle demonstrates that the enclosed area remains finite even though the length becomes infinite as  n -> ∞.

Show[
        GraphicsArray[{DisplayKoch[triangle, 4, 
              DisplayFunction -> Identity], 
            Show[{DisplayKoch[triangle, 4, ShowIterations -> True, 
                  DisplayFunction -> Identity], 
                Graphics[Circle[{0.5, ((Sqrt[3])/2) - (1/(Sqrt[3]))}, (1/(Sqrt[3]))]]}, 
              DisplayFunction -> Identity]}]];

The area is computed by adding or subtracting the area of each new triangle to or from the original area, such that A = (Subscript[A, 0]) ± δ A where the sign is positive for outward or negative for inward triangulation.  The area of an equilateral triangle with sides of length (Subscript[ℓ, n]) is simply ((Sqrt[3])/4) ((Subscript[ℓ, n])^2) so that change in area for iteration becomes δ(Subscript[A, n]) = 
        (Subscript[n, 0]) ((Sqrt[3])/4) (∑_(n = 1)^∞) 
                  (4^(n - 1)) ((Subscript[ℓ, n])^2).  Recognizing that \(TraditionalForm\`ℓ\_\(\(n\)\(\ \)\) = \(\((1\/3)\)\^n\) ℓ\_0\), we find that \(TraditionalForm\`δ A\_n = \(n\_0\) \(ℓ\_0\^2\) \(\@3\/16\) \((4\/9)\)\^n\) is a simple geometric series.

δA = (Subscript[n, 0]) ((Subscript[ℓ, 0])^2) ((Sqrt[3])/16) (∑_(n = 1)^∞)((((4/9)))^n)

For example, the area enclosed by a Koch snowflake obtained by outward triangulation of an equilateral triangle is larger by a factor

((Subscript[A, 0]) + δA)/
        (Subscript[A, 0]) /. {(Subscript[n, 0]) -> 3, (Subscript[ℓ, 0]) -> 1, 
        (Subscript[A, 0]) -> ((Sqrt[3])/4)}

compared with its original area.  It is instructive to compare this area with that of the circumscribed circle of radius (Subscript[ℓ, 0])/(Sqrt[3])

((Subscript[A, 0]) + δA)/(π (r^2)) /. {(Subscript[n, 0]) -> 
            3, (Subscript[ℓ, 0]) -> 1, (Subscript[A, 0]) -> ((Sqrt[3])/4), 
          r -> 1/(Sqrt[3])}  // N

Hence, the snowflake occupies 66% of the circumscribed circle.  It is easy to show that as (Subscript[n, 0]) increases, the area enclosed by the outward Koch curve approaches that of the circumscribed circle as n -> ∞.

Fractal dimension

The size of an object can be measured in terms of an elementary volume (ε^d), where ε is a length scale and d is its dimensionality.  The d-dimensional volume is then the number of elementary volumes, N[ε], contained within the object.  Ordinarily one expects that upon reduction of the length scale, \(TraditionalForm\`ε -> λ\ ε\), the number of elementary volumes scales according to \(TraditionalForm\`N ∝ ε\^\(-d\) ∝ λ\^\(-d\)\).  For example, as the number of sides increases a regular polygon becomes an increasingly accurate approximation to a circle.  If we decrease the length of each side by a rational factor λ and correspondingly increase the number of sides by the factor (λ^(-1)), the length of the perimeter remains constant as the area approaches that of the circumscribed circle.  Thus, a circle is a one-dimensional object of finite length enclosing finite area.  By contrast, the Koch curve also encloses a finite area but has infinite length and, hence, cannot be described as a one-dimensional object.  Therefore, it is necessary to generalize our notions of dimension to include the possibility of nonintegral or fractal dimension.  The procedure proposed by Hausdorff and Besicovitch applies essentially the same scaling argument to objects generated by iterated transformations.

The number of line segments in a Koch curve scales as N ∝ (4^n) while the elementary length scales as λ = ((((1/3)))^n).  Thus, equating these scaling relationships and solving for the dimensionality d

Solve[(4^n) == ((((1/3)))^(-n d)), d] // 
      PowerExpand

% // N

we find that the Koch curve is characterized by a non-integer dimension, d ≈ 1.26.  This fractal dimension occurs because the number of segments increases more rapidly than the length of each segment decreases.  Other transformations which share this behavior will also yield fractal dimensions 1 < d < 2 for which a finite area is enclosed by an infinite length.  Thus, it is useful to distinguish between topological dimension, which is unity for a curve constructed by applying a linear self-similarity transformation to a line, and Hausdorff dimension.  These concepts can be readily generalized to figures constructed by applying self-similarity transformations to simple objects of larger, but integral, topological dimension.  When the Hausdorff dimension d of a closed figure is greater than its topological dimension (Subscript[d, 0]), an infinite (Subscript[d, 0])-dimensional length encloses a finite ((Subscript[d, 0]) + 1)-dimensional area.  Conversely, when d < (Subscript[d, 0]), a finite (Subscript[d, 0])-dimensional length is reduced to a ((Subscript[d, 0]) - 1)-fractal dust (see Cantor set problem below).

This property of fractal curves is sometimes described by analogy with coastlines.  Using a coarse length scale, such as kilometers, the length of a coast has one apparent value.  However, upon closer examination, use of a meter stick would result in a greater length estimate as irregular features of smaller scale are probed.  Decreasing the length scale to centimeter or millimeter scales produces ever greater estimates of the length as finer structures become apparent.  As this process is repeated, the coast line appears to have almost infinite length by the time we reach atomic or subatomic dimensions.  Evidently, coast lines are fractal by nature.