CHAPTER 10, Fractals

10.1 Fractals

Many objects in nature have been modelled mathematically by approximating the curves or surfaces involved by 'smooth' curves or surfaces, in order that the calculus be applied to their study. However, in recent years it has been realized that for certain of these objects, the calculus is not the best tool for their study.
     For example, the motion of a particle suspended in a fluid (Brownian motion), the length of the coastline of an island or the surface area of the human lung. The length of a coast line is dependent on how carefully it is measured; for such an object the more  closely you look at it, the more irregular it appears and the greater the length appears to be. The study of such objects has resulted in a new area of mathematics called Fractal Geometry. We shall use Mathematica to study some of these fractals.

Clear[initialize]
initialize[start_List:{0.`,0.`},dim_Integer:1000]:=(path=Table[Null,{dim}];X=N[start];U={1.`,0.`};path[[1]]=X;count=1)
right[a_]:=U={{t1=Cos[aa=N[a °]],t2=Sin[aa]},{-t2,t1}}.U
left[a_]:=right[-a]
forward[s_]:=path[[++count]]=X+=s U
back[s_]:=forward[-s]
finished:=path=Take[path,count]
showturtlepath:=(finished;Show[Graphics[Line[path]],PlotRange->All,AspectRatio->1])

flakeside[length_,0]:=forward[length]
flakeside[length_,depth_]:=(flakeside[(length/3),depth-1];left[60];flakeside[(length/3),depth-1];right[120];flakeside[(length/3),depth-1];left[60];flakeside[(length/3),depth-1])
KochSnowflake[n_]:=(initialize[];Do[flakeside[1.,n];right[120],{3}] showturtlepath)

10.2 The Koch Snowflake

The Koch Snowflake curve is constructed inductively as follows.
(i) Start with an equilateral triangle with side length = 1, the perimeter of the resulting curve is  L(0) = 3.
   We can use Mathematica to sketch this curve by first activating the above code, and  putting  n = 0, and then activating the command:

KochSnowflake[0]

(ii) Now add equilateral triangles of side length1/3. The resulting curve can be seen by activating the same command but with  n = 1. The perimeter is now  L(1) = 3 * 4* 1/3 = 4.

KochSnowflake[1]

(iii) Add new equilateral triangles as seen by taking  n = 2  in the same code.
The perimeter is now  L(2) = 3 * (4^2) * (1/(3^2)) = (4^2)/3.

KochSnowflake[2]

(iv) Continue in this way, at the n th stage the perimeter is  L(n) = (4^n)/(3^((n-1))). Take in turn  n = 3,4   to get a good idea of what the limiting curve, (the Koch snowflake) looks like. Caution: Mathematica cannot handle values of  n larger than 4 before running out of memory and the program crashing.

KochSnowflake[3]

(v) As we take the limit as  n  tends to infinity, it can be shown that the resulting curve exists, and that since  L(n)  tends to infinity as  n  tends to infinity, this curve has infinite length. Note also that since it is contained in a circle of radius 1, it encloses a finite area.

10.3 Fractal Dimension

The idea of (topological) dimension is intuitively clear. For example, a point or a number of points has 0-dimensions, a curve is 1-dimensional, a surface 2-dimensional, etc. Thus the topological dimension of the snowflake curve is 1.

It turns out that for sets such as the snowflake curve and the Cantor set, there is another very useful idea of dimension, originally called Hausdorff-Besicovitch dimension, now often called Fractal dimension (following B. Mandelbrot). We motivate its definition as follows:

(i)   Given a piece of string, two copies of it result
      in a string twice the size.
(ii)  For a square we need 4 copies.
(iii) For a cube we need 8 copies.
(iv)  For a 4-dimensional cube we need 16 copies.

i.e. to double the (side length) of a d-dimensional cube we need 
               c = (2^d)  copies

so that        d = (Log c)/(Log 2),

Thus for example, if we have an object whose side length doubles if three copies are stuck together, then it would have fractal dimension  ((Log[3])/(Log[2])).

More generally, to increase the d-dimensional cube (object) in size by a factor of  a, we need  c = (a^d)   copies, so

d = ((Log[c])/(Log[d])).

Returning to the snowflake curve, notice that (one side of it) is made up of 4 copies of itself; each 1/3 of the size, so  a = 3, c = 4 and 
                           d = ((Log[4])/(Log[3])) = 1.2616 … 
(Double click on the following command to verify the above assertion).

KochSnowflake[4]

This number  d  is defined to be the fractal dimension of the curve. Roughly speaking a fractalis defined to be a geometrical object whose fractal dimension is strictly greater than its (topological) dimension. The coastline of Britain has  d = 1.25 (approximately), so is very rough. That of the U.S. is closer to 1 (fairly smooth).
   Generally, fractals have the property of being self-similar, i.e. they do not change their appearance significantly when viewed under a microscope of arbitrary magnifying power. This self similarity can take a linear form as in the snow flake curve, or a non linear form as in the Mandelbrott set (see section 6).

10.4 A Fractal Gasket.

Gasket[{a_,b_,c_,d_}]:=With[{ab=(2a+b)/3,ba=(a+2b)/3,bd=(2b+d)/3,db=(b+2d)/3,cd=(2c+d)/3,dc=(c+2d)/3,ac=(2a+c)/3,ca=(a+2c)/3,ad=(2a+d)/3,da=(a+2d)/3,cb=(2c+b)/3,bc=(c+2b)/3},{{a,ab,ac,ad},
{ba,b,bc,bd},{ac,bd,ca,db},{da,db,dc,d},{ca,cb,c,cd}}];
pts={{0,0},{0,1},{1,1},{1,0}}//N;
n=3;
k=Nest[Join @@ Gasket /@ # &, {pts},n]; // AbsoluteTiming
Graphics[{EdgeForm@Black, Polygon@k}]

10.5 A Square Gasket.

startsquare={{0,0},{0,1},{1,1},{1,0}}

{{0,0},{0,1},{1,1},{1,0}}

Show[Graphics[Polygon[startsquare]]]

Gasket[d_]:=Module[{d1,d2,d3,d4,d5},
d1={d[[1]],d[[2]]/3,d[[3]]/3,d[[4]]/3};
d2={d1[[1]]+{0,2/3},d1[[2]]+{0,2/3},d1[[3]]+{0,2/3},d1[[4]]+{0,2/3}};
d3={d1[[1]]+{2/3,2/3},d1[[2]]+{2/3,2/3},d1[[3]]+{2/3,2/3},d1[[4]]+{2/3,2/3}};
d4={d1[[1]]+{1/3,1/3},d1[[2]]+{1/3,1/3},d1[[3]]+{1/3,1/3},d1[[4]]+{1/3,1/3}};
d5={d1[[1]]+{2/3,0},d1[[2]]+{2/3,0},d1[[3]]+{2/3,0},d1[[4]]+{2/3,0}};
]

squares=Gasket[startsquare]

Show[Graphics[Polygon/@squares]]

10.6 The Sierpinski Carpet.

SierpinskiCarpet[{a_,b_,c_,d_}]:=With[{ab=(2a+b)/3,ba=(a+2b)/3,bd=(2b+d)/3,db=(b+2d)/3,cd=(2c+d)/3,dc=(c+2d)/3,ac=(2a+c)/3,ca=(a+2c)/3,ad=(2a+d)/3,da=(a+2d)/3,cb=(2c+b)/3,bc=(c+2b)/3},{{a,ab,ac,ad},{ab,ba,bd,ac},
{ba,b,bc,bd},{bd,bc,cb,ca},{ad,ac,db,da},{da,db,dc,d},{db,ca,cd,dc},{ca,cb,c,cd}}];
pts={{0,0},{0,1},{1,1},{1,0}}//N;
n=3;
k=Nest[Join @@ SierpinskiCarpet /@ # &, {pts},n]; // AbsoluteTiming
Graphics[{EdgeForm@Black, Polygon@k}]

10.7 The Sierpinski Triangle.

Here are two different ways for producing a version of the Sierpinski triangle.

Sierpinski[{a_,b_,c_}]:=With[{ab=(a+b)/2,bc=(b+c)/2,ca=(a+c)/2},{{a,ab,ca},{ab,b,bc},{ca,bc,c}}];
pts={{0,0},{1,0},{1/2,Sqrt[3]/2}}//N;
n=8;
d=Nest[Join @@ Sierpinski /@ # &, {pts},n]; // AbsoluteTiming
Graphics[{EdgeForm@Black, Polygon@d}]

ArrayPlot/@NestList[ArrayFlatten[{{#,0},{#,#}}]&,{{1}},4]

ArrayPlot[Nest[ArrayFlatten[{{#,0},{#,#}}]&,{{1}},8]]

Here is some simple code for producing a version of the Sierpinski carpet.

ArrayPlot/@NestList[ArrayFlatten[{{#,#,#},{#,0,#},{#,#,#}}]&,{{1}},4]

ArrayPlot[Nest[ArrayFlatten[{{#,#,#},{#,0,#},{#,#,#}}]&,{{1}},6]]

This type of construction may be varied to give a wide variety of fractals.

ArrayPlot/@NestList[ArrayFlatten[{{#,0,#},{0,#,#},{#,0,#}}]&,{{2}},6]

10.8 Space Filling Curves

Related to the idea of fractal is the notion of space filling curve. Such curves are sometimes called Peano curves after Giuseppe Peano, who in 1890 discovered the amazing fact that one could map an interval continuously onto a square. By definition,a curve, say in Euclidean space is a continuous function defined on an interval with range in the space. In practice we usually identify the range of the function with the curve. The Koch Snowflake curve can be thought of as a curve in this sense. Peano defined his space filling curve as the limiting curve in a certain sequence of curves, and then showed that the limit curve had the required properties. We use Mathematica to construct the first five iterations of Hilbert's version of a space filling curve. First activate the code below and then double click on Hilbert[n], taking successively the values for n, n=1,2,3,4,5. (If you exceed n=5 the program may crash).

BeginPackage["RecursiveTurtle`"];
RecursiveTurtle::"usage"="RecursiveTurtle[recursion, axiom, depth, angle, steplength, startposn:{0.,0.}, startdir:{1.,0.}, dim:1000] returns the image obtained by applying a turtle to the sequence of movements generated in depth many rewritings of the axiom. Characters are F, B, -, +. If the path is expected to  visit more than 1000 points, a larger value of dim should be used. The rewrite rules in recursion should each transform a  single character."
path::"usage"="path contains the turtle's coordinates  and may be used externally";

"RecursiveTurtle[recursion, axiom, depth, angle, steplength, startposn:{0.,0.}, startdir:{1.,0.}, dim:1000] returns the image obtained by applying a turtle to the sequence of movements generated in depth many rewritings of the axiom. Characters are F, B, -, +. If the path is expected to  visit more than 1000 points, a larger value of dim should be used. The rewrite rules in recursion should each transform a  single character."

Begin["RecursiveTurtle`Private`"];
RecursiveTurtle[recursion_List,axiom_String,depth_Integer,angle_,steplength_,startposn_List:{0,0},startdir_List:{1,0},dim_:1000]:=(symboltable={"F"->forward,"B"->back,"-"->right,"+"->left};X=N[startposn];U=N[startdir];count=1;c=Cos[N[angle Degree]];s=Sin[N[angle Degree]];Attributes[turtle]=Listable;turtle[forward]:=path[[++count]]=X+=steplength U;turtle[back]:=path[[++count]]=X-=steplength U;turtle[left]:=U=rotateleft.U;turtle[right]:=U=rotateright.U;replacerules=(First[Characters[#1[[1]]]]->Characters[#1[[2]]]&)/@recursion/.symboltable;path=Table[Null,{dim}];path[[1]]=X;rotateleft={{c,-s},{s,c}};rotateright=Transpose[rotateleft];turtle[Nest[#1/.replacerules&,Characters[axiom]/.symboltable,depth]];path=Take[path,count];Show[Graphics[Line[path]],AspectRatio->1,Axes->None,PlotRange->All])
End[];EndPackage[];

Hilbert[n_]:=Show[RecursiveTurtle[{"X"->"-YF+XFX+FY-","Y"->"+XF-YFY-FX+"},"Y",n,90,(2.^(-n)),(4^n)],AspectRatio->1,PlotRange->{{0,1},{0,1}}]

Hilbert[5]