Self-Similar Fractals
by Brian H. Marston
http://www.fatdays.com/
webguy@fatdays.com

Initializations

Off[General::spell];
Needs["Graphics`Colors`"]
Needs["Graphics`Polyhedra`"]
Off[ParametricPlot3D::ppcom];
Off[ParametricPlot::ppcom];
Off[Plot::plnr];
SetOptions[Limit,Analytic->True];

Introduction

"Fractals are about looking closely and seeing more.  Fractals have to do with bumps that have bumps, cracks that have crookednesses within crookednesses, and atoms that turn out to be universes.  Fractals have to do with the rich structure of our universe that spans all scales from the uncountable galaxies at unthinkable distances to the mysterious inner electric flashes and vibrations of the subatomic realm."
	--Tim Wegner and Bert Tyler, Fractal Creations,  2nd ed., p. 13

The Modeling Power of Fractals

Many objects in nature are so complicated and irregular that they cannot be modeled well using conic sections, polygons, spheres and the other familiar objects of classical geometry.  For example, circulatory systems, clouds, trees, mountains, and coastlines cannot be reduced to combinations of simple shapes from classical geometry.  Where classical geometry ends as a tool for analyzing the complexity of natural objects, fractal geometry begins.  Today, fractals are used to model a wide range of biological and topographical entities and to produce ultra-realistic special effects for movies and video games.

How Long is the Coastline of Britain?

The question "How long is the coastline of Britain?" posed by Benoit Mandelbrot, the father of modern fractal theory, in his book The Fractal Geometry of Nature is not as simple as it appears.  The problem is that one's answer to this question depends on the length of the ruler one uses.  Unlike circles and the other shapes from classical geometry, coastlines are very irregular.  They're full of inlets, bays, and rocky shores.  A shorter measuring stick will fit more snugly in these nooks and crannies and increase the estimated length of the coastline.  Hence, if we measure the length of Britain's coastline using a mile-long ruler, we will get one value.  If we use a shorter ruler, say a yardstick, we will get a larger value because a yardstick can more closely approximate Britain's convoluted boundary.  In fact, as the scale of measurement decreases, the estimated length increases without limit.  Thus, as the length of the ruler approaches zero, the estimated length of the coastline approaches infinity.  This difficulty in measuring due to the irregularity of the object being measured is characteristic of fractal curves and surfaces.

Self-Similarity

Fractal dimension

Calculating Fractal Dimension

Clear[points,line,n]
points=ListPlot[{{0,0},{1/3,0},{2/3,0},{1,0}},
	PlotStyle->{Red,PointSize[.02]},Axes->False,
	DisplayFunction->Identity];
line=Graphics[Line[{{0,0},{1,0}}]];
Show[{line,
	   points,
	   Graphics[Table[
	   			  Text["1/3",{n,.01}],{n,1/6,5/6,1/3}]]},
	   PlotLabel->FontForm["N=3, r=1/3, Nr=1",
	 		{"Times-Bold",14}],
	 	PlotRange->{-.1,.1},
	 	DisplayFunction->$DisplayFunction]

Clear[i,j,n,horizontal,vertical]
horizontal=Table[Graphics[Line[{{0,j},{1,j}}]],
	{j,0,1,1/3}];
vertical=Table[Graphics[Line[{{i,0},{i,1}}]],
	{i,0,1,1/3}];
Show[{
	  horizontal,
	  vertical,
	  Graphics[Table[
	  		Text["1/3",{-.05,n}],{n,1/6,5/6,1/3}]],
	  Graphics[Table[
	  		Text["1/3",{n,-.05}],{n,1/6,5/6,1/3}]]},
	  AspectRatio->1,
	  PlotLabel->FontForm["N=9, r=1/3, Nr^2=1",
	 		{"Times-Bold",14}]]

Clear[i,j,k,n]
Show[{
	Table[Graphics3D[Cuboid[{i,j,k},
	{i+1/3,j+1/3,k+1/3}]],
	{i,0,2/3,1/3},{j,0,2/3,1/3},{k,0,2/3,1/3}],
	Graphics3D[Table[
		Text["1/3",{n,-.05,-.05}],{n,1/6,5/6,1/3}]],
	Graphics3D[Table[
		Text["1/3",{1.05,n,-.05}],{n,1/6,5/6,1/3}]],	
	Graphics3D[Table[
		Text["1/3",{-.05,-.05,n}],{n,1/6,5/6,1/3}]]},
	Boxed->False,
	PlotLabel->FontForm["N=27, r=1/3, Nr^3=1",
	 	{"Times-Bold",14}]]

Solve[n r^d==1,d]

Clear[dimension,d]
dimension[n_Integer?Positive,r_?Positive]:=
	Module[{d},
		d=Log[1/n]/Log[r]//N]

Koch Snowflake

Construction

(*  Given a non-negative integer n, the function 
    Snowflake generates the nth level approximation 
    to the Koch Snowflake. Code from the Notices of 
    the AMS Vol.39 #7 Sept 92, page 709              *)
    
Clear[n,Snowflake,start,finish,doline]
Snowflake[n_Integer?NonNegative]:=
  Show[Graphics[
    Nest[ (#1/.Line[{start_,finish_}] :>doline[start,
             finish]) &,
          {Line[{{0,0},{1/2,Sqrt[3]/2}}],
           Line[{{1/2,Sqrt[3]/2},{1,0}}],
           Line[{{1,0},{0,0}}]},
          n]],
   AspectRatio->Automatic,PlotRange->All]

doline[start_,finish_]:=
  Module[{vec,normal},
  vec=finish-start;
  normal=Reverse[vec] {-1,1} Sqrt[3]/6;
  {Line[{start,start + vec/3}],
  Line[{start + vec/3,start + vec/2 + normal}],
  Line[{start + vec/2 + normal, start + 2 vec/3}],
  Line[{start + 2 vec/3, finish}]
  }
  ];

Snowflake[0]

Snowflake[1]

Snowflake[2]

Snowflake[3]

Snowflake[4]

Snowflake[5]

Dimension

dimension[4,1/3]

Boundary length

Limit[(4/3)^n,n->Infinity]

Sierpinski Gasket

Construction

(* Given a level, the function Gasket generates all the 
   levels up to and including that level in the construction
   of the Sierpinski Gasket, starting with level 0.
   Code adapted from the code written by Patrick 
   Haggerty to produce the graphics for the Sierpinski
   Gasket in Dr. Richard Crownover's book Fractals and
   Chaos.                                                   *)
  
Clear[Gasket,disp,newF,G,F,sierp,i,j,k,w,trans,group1,group2];
Gasket[level_Integer?NonNegative]:=
	Module[{newF,G,F,sierp,i,j,k,trans,group1,group2},
	
	(* trans={{a1,b1,c1,d1,e1,f1},
	   		   {a2,b2,c2,d2,e2,f2},
	   		   {a3,b3,c3,d3,e3,f3}} represents the
	   coefficients of the affine transformations of
	   the iterated function system given by
	   Ti(x)=|ai bi|x + |ei|   i=1,2,3 
 	         |ci di|    |fi|
 	   which produces the Sierpinski Gasket          *) 
	trans={{1/2,0,0,1/2,0,0},{1/2,0,0,1/2,1/2,0},
		{1/2,0,0,1/2,1/4,N[Sqrt[3]]/4}};
	
	group1[list_]:={{list[[1]],list[[2]]},{list[[3]],
		list[[4]]}};
	group2[list_]:={{list[[5]]},{list[[6]]}};
	
	disp[0]=Show[
  		Graphics[Polygon[{{0,0},{1/2,Sqrt[3]/2},{1,0}}]],
  		AspectRatio->Automatic];
	
	F={{0,0},{1/2,Sqrt[3.]/2},{1,0},{0,0}};
	G=F;

	For[k=1,k<level+1,k++,
		F=Partition[Flatten[G],2];
		G={};
		For[i=1,i<(Length[trans]+1),i++,
			newF=F;
			For[j=1,j<(Length[F]+1),j++,
			  newF=ReplacePart[newF,
			    group1[trans[[i]]].F[[j]]+
			    group2[trans[[i]]],j];];
		G=Append[G,newF];];
	G=Partition[Partition[Flatten[G],2],4];
	sierp=Table[Graphics[Polygon[G[[n]]]],{n,1,Length[G]}];
	disp[k]=Show[sierp,AspectRatio->Automatic];];
	;]

Gasket[5]

Show[GraphicsArray[{Table[disp[n],{n,0,2}],
	Table[disp[n],{n,3,5}]}]];

Dimension

dimension[3,1/2]

Area Measure

Clear[n,a,k]
a[k_]:=3^(k-1)/(4)^k;1/(1-x);  
s[n_]:=Sum[a[k],{k,1,n}]//N;
Table[{n,s[n]},{n,10,80,10}]//TableForm

Menger Sponge

Construction

(* Given a level as input, the function Sponge outputs a
picture of the Menger Sponge at that level of construction.
Code adapted from code written by Robert M. Dickau available
at MathSource (http://www.wri.com/WWWDocs/mathsource/)      *)

Clear[Sponge,level,iters,side,cubmat,i,j,k,n,faces]
Sponge[level_Integer?NonNegative]:=
	Module[{iters,side,cubmat,i,j,k,n,faces},
	
	iters = level;
	
	side = 3. ^ iters ;
	
	cubmat (* cuboid-matrix *) =

        Table[
        
        If[i==side + 1. || j==side + 1. || k==side + 1.,
                
                (* Pad the table's edges with zeroes; if you want
                   to see the complement of the sponge, transpose
                   the 0. and 1. directly below. *)
                   
                0., 1.],
                
                {i,1.,side + 1.},{j,1.,side+1.},{k,1.,side+1.}];
	Do[ If[ 
        
        (Mod[Round[i/3.^n + 0.5],3]==2 &&
                (Mod[Round[j/3.^n + 0.5],3]==2 ||
                        Mod[Round[k/3.^n + 0.5],3]==2)) ||
                        
        (Mod[Round[j/3.^n + 0.5],3]==2 &&
                (Mod[Round[i/3.^n + 0.5],3]==2 ||
                        Mod[Round[k/3.^n + 0.5],3]==2)) ||
                        
        (Mod[Round[k/3.^n + 0.5],3]==2 &&
                (Mod[Round[i/3.^n + 0.5],3]==2 ||
                        Mod[Round[j/3.^n + 0.5],3]==2)),
         
         (* then--taking advantage of eightfold symmetry--... *)
         
        (cubmat[[i,j,k]]=0.;
         cubmat[[side+1-i,j,k]]=0.;
          cubmat[[i, side+1-j,k]]=0.;
          cubmat[[i,j,side+1-k]]=0.;
         cubmat[[side+1-i, side+1-j,k]]=0.;
          cubmat[[side+1-i,j,side+1-k]]=0.;
          cubmat[[i, side+1-j, side+1-k]]=0.;
         cubmat[[side+1-i,side+1-j,side+1-k]]=0.;)
         
          (* ...no cuboid goes there *)],
         
                {i,(side+1)/2},{j,(side+1)/2},{k,(side+1)/2},
                	{n,0.,iters-1.}];
	faces = {};

              (* Instead of using the Cuboid graphics primitive,
                 we show only the polygons visible from
                 viewpoints in the default octant. *)

	Do[ 
        If[ cubmat[[i,j,k]]==1. && cubmat[[i,j,k+1.]]==0.
        
        		(* That is, if a face belongs at {i,j,k}
        		   and there's nothing hiding it, add the
        		   appropriate polygon to the list. *),
        
                AppendTo[ faces,
                
                (* cuboid tops... *)
                {{i,j,k+1.},{i,j+1.,k+1.},
                        {i+1.,j+1.,k+1.},{i+1.,j,k+1.}}
                        
                        ] ],
                        
        {i,1.,side},{j,1.,side},{k,1.,side} ];
(* Since the figure looks the same regardless of which axis
   is vertical, the polygon-corner list "faces" is computed
   only for the tops of the cuboids, then rotated twice to get
   lists of sides and fronts. *)
   
	faces = Join[ faces (*tops*),
        (*sides *)	Map[ RotateLeft[#,2]&,faces,{2}],
        (*fronts*)	Map[ 
        	RotateLeft[#,1]*{1,-1,1}+{0,side+2,0}&,
        		faces,{2}]
        			];
	Show[Graphics3D[ {EdgeForm[], Map[ Polygon, faces ]}],
                Boxed->False];
]

Sponge[0]

Sponge[1]

Sponge[2]

Sponge[3] (* Warning: This picture took 2 minutes, 51 seconds
			   to generate on an SGI Indy                      *)

Dimension

dimension[20,1/3]

Volume Measure

Clear[n,a,k]
a[k_]:=7(20^(k-1)/27^k)
s[n_]:=Sum[a[k],{k,1,n}]//N;
Table[{n,s[n]},{n,5,60,5}]//TableForm

1/(1-x) = 1 + x + x2 + ...,  if |x|<1.

Problems

1

Clear[level]
level[0]=Show[
	Graphics[Line[{{0,0},{1,0}}]],
	Epilog->{
		{PointSize[.025],Red,Point[{0,0}]},
		{PointSize[.025],Red,Point[{1,0}]},
		{Text["0",{0,-.1}]},	
		{Text["1",{1,-.1}]}
		},	
	DisplayFunction->Identity];
level[1]=Show[
	Graphics[{Line[{{0,0},{1/3,0}}],
		Line[{{2/3,0},{1,0}}]}],
	Epilog->{
		{PointSize[.025],Red,Point[{0,0}]},
		{PointSize[.025],Red,Point[{1/3,0}]},
		{PointSize[.025],Red,Point[{2/3,0}]},
		{PointSize[.025],Red,Point[{1,0}]},
		{Text["0",{0,-.1}]},	
		{Text["1/3",{1/3,-.1}]},
		{Text["2/3",{2/3,-.1}]},	
		{Text["1",{1,-.1}]}
		},	
	DisplayFunction->Identity];
level[2]=Show[
	Graphics[{
		Line[{{0,0},{1/9,0}}],
		Line[{{2/9,0},{1/3,0}}],
		Line[{{2/3,0},{7/9,0}}],
		Line[{{8/9,0},{1,0}}]}],
	Epilog->{
		{PointSize[.025],Red,Point[{0,0}]},
		{PointSize[.025],Red,Point[{1/9,0}]},
		{PointSize[.025],Red,Point[{2/9,0}]},
		{PointSize[.025],Red,Point[{1/3,0}]},
		{PointSize[.025],Red,Point[{2/3,0}]},
		{PointSize[.025],Red,Point[{7/9,0}]},
		{PointSize[.025],Red,Point[{8/9,0}]},
		{PointSize[.025],Red,Point[{1,0}]},
		{Text["0",{0,-.1}]},	
		{Text["1/9",{1/9,.1}]},
		{Text["2/9",{2/9,-.1}]},
		{Text["1/3",{1/3,.1}]},
		{Text["2/3",{2/3,-.1}]},	
		{Text["7/9",{7/9,.1}]},	
		{Text["8/9",{8/9,-.1}]},
		{Text["1",{1,.1}]}
		},	
	DisplayFunction->Identity];
Show[GraphicsArray[{{level[0]},{level[1]},
	{level[2]}}]]

Answer and Discussion

dimension[2,1/3]

2

a)

Clear[level]
level[0]=ListPlot[{
	{0,0},{1,0}},
	PlotJoined->True,
	PlotStyle->{Red,Thickness[.01]},
	Ticks->{{{.25,"1/4"},{.5,"1/2"},{.75,"3/4"},{1,"1"}},
		{{.25,"1/4"},{-.25,"-1/4"}}},
	DisplayFunction->Identity
	]
level[1]=ListPlot[{
	{0,0},{.25,0},{.25,.25},{.5,.25},{.5,0},{.5,-.25},
		{.75,-.25},{.75,0},{1,0}},
	PlotJoined->True,
	PlotStyle->{Red,Thickness[.01]},
	Ticks->{{{.25,"1/4"},{.5,"1/2"},{.75,"3/4"},{1,"1"}},
		{{.25,"1/4"},{-.25,"-1/4"}}},
	AspectRatio->Automatic,
	DisplayFunction->Identity
	]
Show[GraphicsArray[{level[0],level[1]},
	DisplayFunction->$DisplayFunction]]

dimension[8,1/4]

b)

Clear[level]
level[0]=ListPlot[{
	{0,0},{1,0}},
	PlotJoined->True,
	PlotStyle->{Red,Thickness[.01]},
	DisplayFunction->Identity
	]
level[1]=ListPlot[{
	{0,0},{0,1/6},{1/6,1/6},{1/6,1/3},{1/3,1/3},
	{1/2,1/3},{1/2,1/6},{1/3,1/6},{1/3,0},{1/2,0},
	{2/3,0},{2/3,-1/6},{1/2,-1/6},{1/2,-1/3},{2/3,-1/3},
	{5/6,-1/3},{5/6,-1/6},{1,-1/6},{1,0}},
	PlotJoined->True,
	PlotStyle->{Red,Thickness[.01]},
	DisplayFunction->Identity
	]
Show[GraphicsArray[{{level[0],level[1]}},
	DisplayFunction->$DisplayFunction]]

dimension[18,1/6]

3

Clear[level]
level[0]=Graphics[
	Rectangle[{0,0},{1,1}],
	AspectRatio->1,
	Axes->True,
	DisplayFunction->Identity];
level[1]=Graphics[{
	Rectangle[{0,0},{1/3,1}],
	Rectangle[{2/3,0},{1,1}],
	Rectangle[{1/3,2/3},{2/3,1}],
	Rectangle[{1/3,0},{2/3,1/3}]},
	AspectRatio->1,
	Axes->True,
	DisplayFunction->Identity]	
Show[GraphicsArray[{level[0],level[1]}],
	DisplayFunction->$DisplayFunction]

Answer:

dimension[8,1/3]

Other Fractal Resources

Frequently Asked Questions (FAQ) file for the sci.fractals newsgroup

Fractals and Chaos by Dr. Richard M. Crownover

Fractal Creations, 2nd Ed. by Tim Wegner and Bert Tyler