Fractals

Contents

	Home
$Math Artwork$	Math Artwork
	Fractals
	Fluid Motion
	Physics Simulations
	Engineering
	Minimal Surfaces
	Rendered Artwork
	Hand-made Artwork
	Bug Collection
	Programming
	High Voltage
	Physics Experiments
	The Bible
	Links

Fractals are one of the most interesting puzzles of mathematics. They are made by a simple formula, yet they have such beautiful and complex designs. This page provides simple Mathematica code for many of the most interesting types of fractals. Some of this code may not be very efficient, but it should be enough to give you a basic understanding of the mathematics involved.

Kleinian Double Cusp Group - calculated in Mathematica 4.2, rendered in POV-Ray 3.6.1, 8/10/05

This is my attempt to create the double cusp group on page 269 of Indra’s Pearls. Thanks to Dr. William Goldman for helping me get started with this. See also my Double Spiral. Click here to see some POV-Ray code. Some Mathematica code is shown below. You can also see this code on Roger Bagula’s web site:
(* runtime: 0.06 second *)
ta = 1.958591030 - 0.011278560I; tb = 2; tab = (ta tb + Sqrt[ta^2tb^2 - 4(ta^2 + tb^2)])/2; z0 = (tab - 2)tb/(tb tab - 2ta + 2I tab);
b = {{tb - 2I, tb}, {tb, tb + 2I}}/2; B = Inverse[b]; a = {{tab, (tab - 2)/z0}, {(tab + 2)z0, tab}}.B; A = Inverse[a];
Fix[{{a_, b_}, {c_, d_}}] := (a - d - Sqrt[4 b c + (a - d)^2])/(2 c); ToMatrix[{z_, r_}] := (I/r){{z, r^2 - z Conjugate[z]}, {1, -Conjugate[z]}};
MotherCircle[M1_, M2_, M3_] := ToMatrix[{x0 +I y0, r}] /. Solve[Map[(Re[#] - x0)^2 + (Im[#] - y0)^2 == r^2 &, Fix /@ {M1, M2, M3}], {x0, y0, r}][[2]];
C1 = MotherCircle[b, a.b.A, a.b.A.B]; C2 = MotherCircle[b.a.a.a.a.a.a.a.a.a.a.a.a.a.a.a, a.b.a.a.a.a.a.a.a.a.a.a.a.a.a.a, a.b.A.B];
Reflect[C_, M_] := M.C.Inverse[Conjugate[M]];
orbits = Join[Reverse[NestList[Reflect[#, a] &, C1, 63]], Drop[NestList[Reflect[#, A] &, C1, 63], 1], Reverse[NestList[Reflect[#, a] &, C2, 71]], Drop[NestList[Reflect[#, A] &, C2, 56], 1]];
Show[Graphics[MapIndexed[({{a, b}, {c, d}} = #1; {Hue[#2[[1]]/15], Disk[{Re[a/c], Im[a/c]}, Re[I/c]]}) &, orbits]], PlotRange -> 35{{-1, 1}, {-1, 1}}, AspectRatio -> Automatic]

This animation shows the set morphing into a single cusp group.

Links
Explanation - by David Wright, coauthor of Indra’s Pearls
Kleinian Gallery - beautiful Kleinian groups by Jos Leys
Limit Sets - interesting animations by Jeffrey Brock, see his 3D bending
Indra’s Pearls course - Kleinian groups with David Wright
Swirlique - program for drawing Kleinian limit sets
Gallery - includes some nice Kleinian groups, by Curtis McMullen
Double Spiral - POV-Ray rendering by Jeffrey Pettyjohn
Doyle spirals - hexagonal circle packings by Jos Leys
Jellyfish movie - by Curtis McMullen
Kreistiefe - hexagonal circle packings Java applet, by Tim Hoffmann and Nadja Kutz
Hausdorff dimension - extended non-negative real number associated to any metric space

Kleinian Quasifuchsian Limit Set - calculated in Mathematica 4.2, rendered in POV-Ray 3.6.1, 4/15/06

Here is a Sunset Moth “blown about” inside a Quasifuchsian limit set. Originally, Felix Klein described these fractals as “utterly unimaginable”, but today we can visualize these fractals with computers. Here is some Mathematica code:
(* runtime: 12 seconds *)
ta = tb = 1.91 + 0.05I; tab = (ta tb + Sqrt[ta^2tb^2 - 4(ta^2 + tb^2)])/2; z0 = (tab - 2)tb/(tb tab - 2ta + 2I tab);
b = {{tb - 2I, tb}, {tb, tb + 2I}}/2; B = Inverse[b]; a = {{tab, (tab - 2)/z0}, {(tab + 2)z0, tab}}.B; A = Inverse[a];
Reflect[{{a_, b_}, {c_, d_}}, z_] := (b + a z)/(d + c z);
ReflectList[C_, zlist_] := Reflect[C, #] & /@zlist; Children[zlist_] := ReflectList[#, zlist] & /@ {a, b, A, B};
zlist = {0.23 + 0.03 I, 0.18 + 0.05 I, 0.62 + 0.45 I, 0.86 + 0.73 I, 0.91 + 0.89 I, 0.88 + 0.97 I, 0.75 + 0.98 I, 0.48 + 0.88 I, 0.25 + 0.85 I, 0.04 + 0.79 I, -0.02 + 0.67 I, -0.1 + 0.78 I, -0.14 + 0.77 I, -0.24 + 0.84 I, -0.24 + 0.77 I, -0.41 + 0.88 I, -0.39 + 0.77 I, -0.5 + 0.82 I, -0.48 + 0.74 I, -0.82 + I, -0.86 + 0.96 I, -0.68 + 0.79 I, -0.7 + 0.74 I, -0.89 + 0.81 I, -0.74 + 0.64 I, -0.77 + 0.6 I, -0.91 + 0.65 I, -0.8 + 0.51 I, -0.87 + 0.44 I, -0.71 + 0.32 I, -0.44 + 0.19 I, -0.07 + 0.08 I, -0.38 + 0.03 I};
zlists1 = zlists2 = {0.5 + 0.125Join[Reverse[Conjugate[zlist]], zlist]}; test[zlist1_, zlist2_] := Abs[zlist2[[1]] - zlist1[[1]]] < 0.05; While[zlists2 =!= {},zlists2 = Complement[Flatten[Children /@ zlists2, 1], zlists1, SameTest -> test]; zlists1 = Union[zlists2, zlists1, SameTest -> test]];
Show[Graphics[Line[{Re[#], Im[#]} & /@ #] & /@ zlists1], AspectRatio -> Automatic]

Magnetic Pendulum Strange Attractor
Mathematica 4.2 version: 1/24/05; C++ version: 10/27/07; 3D rendered in POV-Ray 3.1

This chaotic strange attractor represent the final resting positions for a magnetic pendulum suspended over some magnets (shown as black dots). It kind of looks like mixed paint. The 2D animation shows what happens as you decrease the damping factor. The 3D animation was shaded by path length. Click here to see an "outrageous" picture of the far field. Here is some Mathematica code:
(* runtime: 25 seconds, increase n for higher resolution *)
n = 40; h = 0.25; g = 0.2; mu = 0.07; zlist = {Sqrt[3] + I, -Sqrt[3] + I, -2I};
image = Table[z2 = z[25] /. NDSolve[{z''[t] == Plus @@ ((zlist - z[t])/(h^2 + Abs[zlist - z[t]]^2)^1.5) - g z[t] - mu z'[t], z[0] == x + I y, z'[0] == 0}, z, {t, 0, 25}, MaxSteps -> 200000][[1]]; r = Abs[z2 - zlist]; i = Position[r, Min[r]][[1, 1]]; Hue[i/3], {y, -5.0, 5.0, 10.0/n}, {x, -5.0,5.0, 10.0/n}];
Show[Graphics[RasterArray[image]], AspectRatio -> 1]

Here is some Mathematica code to numerically solve this using the Beeman integration scheme with the predictor-corrector modification:
(* runtime: 45 seconds, increase n for higher resolution *)
n = 40; tmax = 25; dt = 0.1; h = 0.25; g = 0.2; mu = 0.07; zlist = {Sqrt[3] + I, -Sqrt[3] + I, -2I};
image = Table[z = x + I y; v = a = a1 = 0; Do[z += v dt + (4a - a1)dt^2/6; vpredict = v + (3a - a1)dt/2; a2 = Plus @@ ((zlist - z)/(h^2 + Abs[zlist - z]^2)^1.5) - g z - mu vpredict; v += (2a2 + 5a - a1)dt/6; a1 =a; a = a2, {t, 0, tmax, dt}]; r = Abs[z - zlist]; Hue[Position[r, Min[r]][[1, 1]]/3], {y, -5.0, 5.0, 10.0/n}, {x, -5.0, 5.0, 10.0/n}];
Show[Graphics[RasterArray[image]], AspectRatio -> 1]

The picture on the left shows another version with five magnets. See also my linked pendulum and spherical pendulum.

Links
Magnetic Pendulum - C++ program by Ingo Berg, uses Beeman integration scheme, here is a nice big picture and here is an animation
YouTube Video - explains magnetic pendulum simulation
Butterfly Effect - systems that are highly sensitive to initial conditions
homemade magnetic pendulum - by David Williamson

Quaternion Julia Set Fractal - Mathematica 4.2: 7/9/04, POV-Ray 3.1: 8/5/04

q_n+1 = q_n²+q_c, q_c = -0.2+0.8i+0j+0k
Quaternions are 4 dimensional hypercomplex numbers (as opposed to the usual 2D complex numbers). Here is some Mathematica code demonstrating one way to create a quaternion Julia Set fractal by slowly marching cubes (voxels) toward the boundary:
(* runtime: 2.7 minutes *)
n = 50; qc = {-0.2, 0.8, 0, 0}; Abs2[q_] := Sqrt[Plus @@ Map[#^2 &, q]];
Square[{x_, y_, z_, w_}] := {x^2 - y^2 - z^2 - w^2, 2 x y, 2 x z, 2 x w};
bound[x_, y_, z_] := Module[{q = {x, y, z, 0}, i = 0}, While[i < 12 && Abs2[q] < 2, q = Square[q] + qc; i++]; i == 12];
image = Table[Module[{z = 1.5}, While[z > -1.5 && ! bound[x, y, z], z -= 0.01]; {x, y, z}], {y, -1.5, 1.5, 3/n}, {x, -1.5, 1.5, 3/n}];
ListDensityPlot[Map[#[[3]] &, image, {2}], Mesh -> False, Frame -> False]
We can use fast Phong shading to make it look 3D:
(* runtime: 1 second *)
Normalize[x_] := x/Sqrt[x.x];
ListDensityPlot[Table[Module[{q = image[[i, j]], normal, L, reflect}, normal = Normalize[{(image[[i, j + 1, 3]] - image[[i, j - 1, 3]])q[[3]], (image[[i + 1, j, 3]] - image[[i - 1, j, 3]]) q[[3]], 6/n}]; L = {1, 1, 0} - q; reflect = Normalize[2(normal.L)normal - L]; (normal.Normalize[L] + reflect[[3]])/4 + 1/2], {i, 2, n}, {j, 2, n}], Mesh -> False, Frame -> False]

POV-Ray has an efficient algorithm for ray tracing these fractals:
// runtime: 8 seconds
camera{location <1,1,-2> look_at 0}
light_source{<1.5,1.5,-1.5>,1}
julia_fractal{<-0.2,0.8,0,0> quaternion sqr max_iteration 12 precision 400 pigment{rgb 1}}

Links
Explanation - by Jeff Kelly
Hypercomplex Iterations - nice pictures

Inverse Quaternion Julia Set Fractal
original version: Mathematica 4.2, 7/18/04; animated version: POV-Ray 3.1, 7/4/06

Here is the same fractal again using the inverse Julia set technique. The 4th dimension has been color-coded. This method is much faster, but some regions are faint because they attract much slower. See also my POV-Ray code and rotatable 3D version. Here is some Mathematica code:
(* runtime: 5 seconds *)
Sqrt2[q_] := Module[{r = Sqrt[Plus @@ Map[#^2 &, q]], a, b}, a = Sqrt[(q[[1]] + r)/2]; b = (r - q[[1]]) a/(q[[2]]^2 + q[[3]]^2 + q[[4]]^2); {a, b q[[2]], b q[[3]], b q[[4]]}];
QInverse[qlist_] := Flatten[Map[Module[{q = Sqrt2[# - qc]}, {q, -q}] &, qlist], 1];
qc = {-0.2, 0.8, 0, 0}; qlist = {{1.0, 1.0, 1.0, 1.0}};
Do[qlist = QInverse[qlist], {12}];
Show[Graphics3D[{PointSize[0.005], {Hue[#[[4]]], Point[Delete[#, 4]]} & /@ qlist}, Boxed -> False, Background -> RGBColor[0, 0, 0]]]

Links
QJulia - inverse quaternion Julia set program by Chris Laurel
Quaternion Julia Crystal - laser-etched cube by Bathsheba Grossman

Mandelbrot Set Zoom - original version: Java, 5/24/01; animated version: C++, 2/16/05

z_n+1 = z_n²+z_c
This animation zooms in on the Mandelbrot set by a factor of 10¹⁵. At this high resolution, double precision numbers are inadequate. Therefore, this animation was created using “double-double” precision numbers, adapted from Keith Brigg’s code for double-doubles. See also my Java program and C program. Here is some Mathematica code for this fractal:
(* runtime: 1 minute *)
Mandelbrot[zc_] := Module[{z = 0, i = 0}, While[i < 100 && Abs[z] < 2, z = z^2 + zc; i++]; i];
DensityPlot[Mandelbrot[xc + I yc], {xc, -2, 1}, {yc, -1.5, 1.5}, PlotPoints -> 275, Mesh -> False, Frame -> False, ColorFunction -> (If[# != 1, Hue[#], Hue[0, 0, 0]] &)]

Here is a faster version:
(* runtime: 7 seconds *)
Mandelbrot = Compile[{{zc, _Complex}}, Length[FixedPointList[#^2 + zc &, zc, 100, SameTest -> (Abs[#] > 2 &)]]];
DensityPlot[Mandelbrot[xc + I yc], {xc, -2, 1}, {yc, -1.5, 1.5}, PlotPoints -> 275, Mesh -> False, Frame -> False, ColorFunction -> (If[# != 1, Hue[#], Hue[0, 0, 0]] &)]

POV-Ray has a built-in function for this fractal:
// runtime: 0.5 second
camera{orthographic location <-0.5,0,-1.5> look_at <-0.5,0,0> angle 90}
plane{z,0 pigment{mandel 100 color_map{[0 rgb 0][1/6 rgb <1,0,1>][1/3 rgb 1][1 rgb 1][1 rgb 0]}} finish{ambient 1}}

Links
Deep-zooming - Fractint can zoom in by magnification factors of up to 10¹⁶⁰⁰
Mandelbrot Deep Zoom - zooming animation, (10⁹⁵)

Wada Basin - Mathematica 4.2, 6/19/05; C++ version: 9/11/08

To make a Wada basin, simply put some Christmas tree ornaments together. The following Mathematica code demonstrates how to write a very simple ray tracer. Click here to download a simple C++ program.
(* runtime: 65 seconds, increase n for higher resolution *)
n = 140; r = 1; plight = {0, 0, 1}; Normalize[x_] := x/Sqrt[x.x];
spheres = {{{1, 1.0/Sqrt[3], 0}, {1, 0,0}}, {{-1, 1.0/Sqrt[3], 0}, {0, 1, 0}}, {{0, -2.0/Sqrt[3], 0}, {0, 0, 1}}, {{0, 0, -Sqrt[8.0/3]}, {1, 1, 1}}};
Intersections[p_, dir_] := Module[{a = p.dir, b}, b = r^2 + a^2 - p.p; If[b > 0, Select[Sqrt[b]{-1, 1} - a, # > 10^-10 &], {}]];
RayTrace[p_, dir_, depth_] := Module[{dlist = Map[Min[Intersections[p - #[[1]], dir]] &, spheres], d, i, p2, normal, color, lightdir, intensity}, d = Min[dlist]; If[d < Infinity, p2 = p + d dir; i = Position[dlist, d][[1, 1]]; normal = Normalize[p2 - spheres[[i, 1]]]; color = If[depth > 50, {0, 0, 0}, RayTrace[p2, Normalize[dir - 2(normal.dir)normal], depth + 1]]; lightdir = Normalize[plight - p2]; intensity = normal.lightdir; If[intensity > 0, shadowed = Or @@ Map[Intersections[p2 - #[[1]], lightdir] =!= {} &, spheres]; If[! shadowed,color += intensity spheres[[i, 2]]]]; color, {0, 0, 0}]];
Show[Graphics[RasterArray[Table[RGBColor @@ Map[Min[1, Max[0, #]] &, RayTrace[{0, 0, 1},Normalize[{x,y, -1}], 1]], {y, -0.5, 0.5, 1.0/n}, {x, -0.5, 0.5, 1.0/n}]], AspectRatio -> 1]]

You can also easily ray trace Wada Basins using POV-Ray:
// runtime: 4 seconds
global_settings{max_trace_level 50}
camera{location z look_at 0}
light_source{z,1}
#default{finish{reflection 1 ambient 0}}
sphere{<1,sqrt(1/3),0> 1 pigment{rgb <1,0,0>}}
sphere{<-1,1/sqrt(3),0> 1 pigment{rgb <0,1,0>}}
sphere{<0,-2/sqrt(3),0> 1 pigment{rgb <0,0,1>}}
sphere{<0,0,-sqrt(8/3)> 1 pigment{rgb 1}}

Links
Indra's Deep Sea Life - beautiful POV-Ray images by Kevin Wampler
Scratch a Pixel - ray-tracing tutorial with source code
Javascript Ray Tracer - simple Javascript ray tracer by John Haggerty

Weierstrass Function, POV-Ray 3.6.1: 5/18/08, Mathematica 4.2: 3/2004

The infamous Weierstrass function is an example of a function that is continuous but completely undifferentiable. Here is some Mathematica code:
(* runtime: 0.7 second *)
Plot3D[Sum[Sin[Pi k^a x]/(Pi k^a), {k, 1, 50}], {x, 0, 1}, {a, 2, 3}, PlotPoints -> 100]

Polynomial Roots - Mathematica 4.2, 12/29/08

Strange fractal patterns emerge when you plot the complex roots of high order polynomials. This picture shows all the roots for all possible combinations of 18th order polynomials with coefficients of ±1. Notice that some parts bear a remarkable resemblance to the Heighway Dragon Curve. You can easily find the roots using Mathematica's Root function:
(* runtime: 34 seconds *)
order = 12; n = 275; image = Table[0.0, {n}, {n}];
Do[Do[z = N[Root[Sum[(2Mod[Floor[(t - 1)/2^i], 2] - 1) #^(order - i), {i, 0, order}], root]]; {j,i} = Round[n({Re[z], Im[z]}/1.5 + 1)/2]; If[0 < i <= n && 0 < j <= n, image[[i, j]]++], {root, 1, order}], {t, 1, 2^order}];
ListDensityPlot[image, Mesh -> False, Frame -> False, PlotRange -> {0, 4}]

You can also find all the complex roots of any high order polynomial using the Durand-Kerner method:
(* runtime: 0.1 second *)
n = 18; coefficients = Table[1.0, {n + 1}]; f[z_] := Sum[coefficients[[i + 1]]z^i, {i, 0, n}];
z = Table[(0.4 + 0.9I)^i, {i, 0, n - 1}]; z0 = 0.0z; While[z != z0, z0 = z; z -= f[z]/Table[Product[If[i != j, z[[i]] - z[[j]], 1], {j, 1, n}], {i, 1, n}]]; Chop[z]

Links
Constrained Coefficients - beautiful plots by Loki J�rgenson
description - by John Baez
Integer Coefficients - slightly different plots by Dan Christensen

Nebulabrot - Mathematica 4.2, 8/3/04; C++ version: 5/16/07

This beautiful variation of the Mandelbrot set kind of resembles a nebula or a hologram. The Nebulabrot is a variation of the Buddhabrot technique. Here is some Mathematica code. The coloring method was adapted from Paul Bourke’s web site:
(* runtime: 10 minutes *)
Mandelbrot[zc_] := Length[NestWhileList[#^2 + zc &, 0, Abs[#] < 2 &, 1, 1000]];
n = 275; image = Table[0, {n}, {n}];
Do[If[Mandelbrot[xc + I yc] < 1000, Module[{z = 0, iter = 0}, While[Abs[z] < 2, z = z^2 + xc + I yc; {i, j} = Floor[n{Im[z] + 1.5, Re[z] + 2}/3]; If[0 < i <= n && 0 < j <= n, image[[i, j]] += xc + I yc]; iter++]]], {xc, -2.0, 1.0, 0.01}, {yc, -1.5, 1.5, 0.01}];
Hue2[hue_, x_] := Hue[hue, Min[1, 2(1 - x)], Min[1, 2x]];
Show[Graphics[RasterArray[Map[Hue2[Arg[#]/Pi, Min[Abs[#]/18, 1]] &, image, {2}]], AspectRatio -> 1]]

Diffusion Limited Aggregation (DLA) Fractal - Mathematica 4.2, 8/12/04; C++ version: 2/10/05

This fractal simulates a diffusive growth process similar to that often found in nature. It is generated by single points that randomly drift around until they find something to stick to. The title on this page was generated using this technique. Here is some Mathematica code:
(* runtime: 33 seconds *)
n = 100; x = y = n/2; SeedRandom[0]; image = Table[0, {n}, {n}];
Do[{x, y} = Floor[n(0.5 + 0.1{Cos[theta], Sin[theta]})] + 1; image[[x, y]] = 1, {theta, 0, 2 Pi, Pi/180}];
Do[theta = 2 Pi Random[]; {x, y} = Floor[n(1 + {Cos[theta], Sin[theta]})/2]; drift = True; While[drift, {x, y} = Mod[{x + Random[Integer, {-1, 1}], y + Random[Integer, {-1, 1}]}, n]; drift = Plus @@ Flatten[image[[Mod[x + {-1,0, 1}, n] + 1, Mod[y + {-1, 0, 1}, n] + 1]]] == 0]; image[[x + 1, y + 1]] = 1 - t/360, {t, 0, 360}];
ListDensityPlot[image, Mesh -> False, Frame -> False]

Links
Andy Lomas - impressive 3D aggregation simulations
Lichtenberg Figures - beautiful electric discharge patterns
3D DLA Fractals - by Paul Bourke
Java Applet - by Chi-Hang Lam
Dielectric Breakdown Model - simulates lightning bolts
Snowflake Simulations - rendered in POV-Ray, by Theodore Kim and Ming Lin
Snowfake Movies - cellular automata
Snow Crystals - Cal Tech’s beautiful guide to snowflakes, be sure to see the growing snowflake movie
giant snowflake - from an Alaskan permafrost tunnel
Frost Flowers - very weird ice that looks like hair
rapid crystal growth - as part of the National Ignition Facility Project
Century Egg - unusual Chinese delicacy that can have dendrite-like patterns growing on it

Newton-Raphson Fractal - Mathematica 4.2, 6/17/04

f(z) = z⁵-1, z_n+1 = z_n - f(z_n)/f'(z_n)
It is surprisingly easy to make these beautiful fractals. Here is some Mathematica code:
(* runtime: 25 seconds *)
f[z_] := z^5 - 1;
Show[Graphics[RasterArray[Table[Module[{data = FixedPointList[# - f[#]/f'[#] &, x + I y, 29], z}, z = 1 - Length[data]/30; Hue[Arg[Last[data]]/(2Pi),Min[1, 2(1 - z)], Min[1, 2z]]], {y, -1.5, 1.5, 3.0/275}, {x, -1.5, 1.5, 3.0/275}]]], AspectRatio -> 1]

Links
Newton-Raphson Fractal - Mathematica code by Mark McClure
Newton-Raphson Fractal - Mathematica code by Robert M. Dickau

“The Watering Hole” Orbit Trap - old version: Mathematica 4.2, 2/26/05; animated version: C++, 7/12/06

Orbit trapping is a popular technique for mapping arbitrary shapes to a fractal. This orbit trap is shaped like a Sunset Moth. Here is some basic Mathematica code:
(* runtime: 18 seconds, Note: black areas of Picture.bmp will be transparent *)
image = Import["C:/Picture.bmp"][[1, 1]]/255.0; i2 = Length[image]; j2 = Length[image[[1]]]; n = 275;
map[z_] := Round[{i2(Im[z] + 1), j2(Re[z] + 1)}/2];
trapped[i_, j_] := 0 < i <= i2 && 0 < j <= j2 && image[[i, j]] != {0, 0, 0};
OrbitTrap[zc_] := Module[{z = 0, i, j, k = 0}, While[k < 100 && Abs[z] < 2 && (k < 2 || ! trapped @@ map[z]), z = z^2 + zc;k++]; {i, j} = map[z]; If[trapped[i, j], RGBColor @@ image[[i, j]], RGBColor[0, 0, 0]]];
Show[Graphics[RasterArray[Table[OrbitTrap[xc + I yc], {yc, -1.5, 1.5, 3.0/n}, {xc, -2, 1, 3.0/n}]], AspectRatio -> 1]]

Gears Orbit Trap - POV-Ray 3.1, C++, 7/12/06

Here is the same fractal again using a gear-shaped orbit trap.

Mandelbrot Set Tessellation - Mathematica 4.2, 6/14/04

The left image is distorted because it is not a hyperbolic tiling. For more information, see Bill Rood’s web site and the book “The Colours of Infinity”. Here is some basic Mathematica code:
(* runtime: 70 seconds *)
image = Import["C:/Picture.jpg"][[1, 1]]; n = Length[image]; m = Length[image[[1]]];
Show[Graphics[RasterArray[Table[Module[{zc = x + I y, z = 0.0}, Do[z = z^2 + zc, {25}]; RGBColor @@ (image[[Floor[n Mod[Re[Log[Log[z]]]m/n, 1]] + 1, Floor[m Mod[2Arg[zc], 1]] + 1]]/255.0)], {y, -1.5, 1.5, 3/275}, {x, -2, 1, 3.0/275}]], AspectRatio -> 1]]

Golden Ratio Spiral Orbit Trap - Mathematica 4.2, 3/3/05

Here is an orbit trap in the shape of the Golden Ratio Spiral. Here is some Mathematica code:
(* runtime: 1 minute *)
phi = 0.5(1 + Sqrt[5]);
f[z_] := Module[{r = Log[Abs[z]]/(4Log[phi]) - Arg[z]/(2Pi)}, 18Abs[r - Round[r]]];
OrbitTrap[zc_] := Module[{z = 0, i =0}, While[i < 100 && Abs[z] < 100 && (i < 2 || f[z] > 1), z = z^2 + zc; i++]; Hue[Arg[z]/(2Pi), Min[1, 2 f[z]], Max[0, Min[1, 2(1 - f[z])]]]];
Show[Graphics[RasterArray[Table[OrbitTrap[xc + I yc], {yc, -1.5, 1.5, 3.0/275}, {xc, -2, 1, 3.0/275}]], AspectRatio -> 1]]

Frequency Filtered Random Noise - Mathematica 4.2, 9/5/04

This technique uses a Fourier transform to generate random periodic textures. See also my image deconvolution. Here is some Mathematica code:
(* runtime: 20 seconds *)
n = 275; SeedRandom[1];
ListDensityPlot[Abs[InverseFourier[Fourier[Table[Random[], {n}, {n}]] Table[Exp[-((j/n - 0.5)^2 + (i/n - 0.5)^2)/0.025^2], {i, 1, n}, {j, 1, n}]]], Mesh -> False, Frame -> False, ColorFunction -> (Hue[# + 0.5] &)]

Here is some Mathematica code for an animated version. This technique can be used for simulating water surfaces:
(* runtime: 18 seconds *)
n = 64; SeedRandom[0];
Scan[ListDensityPlot[#, Mesh -> False, Frame -> False, ColorFunction -> (Hue[# + 0.5] &)] &, Abs[InverseFourier[Fourier[Table[Random[], {32}, {n}, {n}]]Table[Exp[-((i/n - 0.5)^2 + (j/n - 0.5)^2 + (t/32 - 0.5)^2)/0.025^2], {t, 1, 32}, {i, 1, n}, {j, 1, n}]]]]

Perlin Noise - Mathematica 4.2, 8/16/04

This is a popular technique for generating all kinds of random periodic textures such as artificial terrain, clouds, water, wood, fire, marble, etc.. It was invented by Ken Perlin. Most random number generators base subsequent random numbers on their previous numbers. In this way, highly complicated “random” scenes can be easily reconstructed from a single initial seed number. Here is some Mathematica code:
(* runtime: 2 minutes *)
Interpolate[{y0_, y1_, y2_, y3_}, mu_] := (y3 - y2 - y0 + y1)mu^3 + (2 y0 - 2 y1 + y2 - y3) mu^2 + (y2 - y0) mu + y1;
m = 14; SeedRandom[0]; noise = Table[Random[], {m}, {m}];
Noise[x0_, y0_] := Module[{x = Mod[x0, 1], y = Mod[y0, 1], i0, j0}, i0 = Floor[m y]; j0 = Floor[m x]; Interpolate[Table[Interpolate[Table[noise[[Mod[i, m] + 1, Mod[j, m] + 1]], {j, j0 - 1, j0 + 2}], m x - j0], {i, i0 - 1, i0 + 2}], m y - i0]];
a = 2; b = 2; n = 275; image = Table[Sum[Noise[b^i x, b^i y]/a^i, {i, 0, 4}], {x, 0.0, 1.0, 1.0/n}, {y, 0.0, 1.0, 1.0/n}];
ListDensityPlot[image, Mesh -> False, Frame -> False]

Here is some Mathematica code that uses a triangle wave to make Perlin noise that looks like marble:
(* runtime: 4 seconds *)
ListDensityPlot[Table[t = 5j/n +image[[i, j]]; Mod[t, 1] - 2Mod[t + 0.5, 1]Floor[Mod[t, 1] + 0.5], {i, 1, n}, {j, 1, n}], Mesh -> False, Frame -> False, ColorFunction -> (RGBColor[1 - 0.5#, 0.5(1 - #), 0] &)]

Artificial Terrain - Mathematica 4.2, 10/17/04

Here is some Mathematica code to plot the Perlin noise as an artificial terrain (you must run the above Perlin noise code before you can run this code):
(* runtime: 10 seconds *)
Gradient2[x_, grad_] := Module[{i = 1, n = Length[grad]}, While[i <= n && grad[[i, 1]] < x, i++]; RGBColor @@ If[1 < i <= n, Module[{x1 = grad[[i - 1, 1]], x2 = grad[[i, 1]]}, ((x2 - x) grad[[i - 1, 2]] + (x - x1)grad[[i, 2]])/(x2 - x1)], grad[[Min[i, n], 2]]]];
EarthTones = {{0, {0, 0.1, 0.24}}, {0.1, {0, 0.6, 0.6}}, {0.1, {0.9, 0.8, 0.6}}, {0.2, {0.5, 0.4, 0.3}}, {0.5, {0.2, 0.3, 0.1}}, {0.8, {0.7, 0.6, 0.4}}, {0.8, {1, 1, 1}}, {1, {1, 1, 1}}};
Show[SurfaceGraphics[Map[Max[#, 1] &, image, {2}], Mesh -> False, Boxed -> False, BoxRatios -> {1, 1, 1/12}, ColorFunction -> (Gradient2[1.3#, EarthTones] &), Background -> RGBColor[0, 0, 0]]]

See also my Terragen terrain.

Spherical Perlin Noise - Mathematica 4.2, 11/9/04

Here is some Mathematica code to map Perlin noise to a sphere as described on Paul Bourke’s web site. This is accomplished by evaluating a 3D Perlin noise “cloud” at points on the sphere:
(* runtime: 7 minutes *)
m = 14; SeedRandom[0]; noise = Table[Random[], {m}, {m}, {m}];
Noise[x0_, y0_, z0_] := Module[{x = Mod[x0, 1], y = Mod[y0, 1], z = Mod[z0, 1], i0, j0, k0}, i0 = Floor[m y]; j0 = Floor[m x]; k0 = Floor[m z]; Interpolate[Table[Interpolate[Table[Interpolate[Table[noise[[Mod[i, m] + 1, Mod[j, m] + 1, Mod[k, m] + 1]], {j, j0 - 1, j0 + 2}], m x - j0], {i, i0 - 1, i0 + 2}], m y - i0], {k, k0 - 1, k0 + 2}], m z - k0]];
a = 2; b = 2; Perlin[x_, y_, z_] := Sum[Noise[b^i x, b^i y, b^i z]/a^i, {i, 0, 4}];
image = Table[Perlin[Cos[theta]Sin[phi], Sin[theta]Sin[phi], Cos[phi]], {phi, 0, Pi, Pi/180}, {theta, 0, 2Pi, Pi/180}];
ListDensityPlot[image, Mesh -> False, Frame -> False, AspectRatio -> Automatic, ImageSize -> {360, 180}]

Artificial Planet - Mathematica 4.2, 11/9/04

Here is some Mathematica code to plot Perlin noise as an artificial planet (you must run the above Perlin noise code before you can run this code):
(* runtime: 7 minutes *)
ParametricPlot3D[Module[{x = Cos[theta]Sin[phi], y = Sin[theta]Sin[phi], z = Cos[phi], w}, w = Max[0, Perlin[x, y, z] - 0.8]; Append[(1 + 0.1w){x, y, z}, {EdgeForm[], Gradient2[w, EarthTones]}]], {theta, 0, Pi}, {phi, 0, Pi}, PlotPoints -> {361, 181}, Background -> RGBColor[0, 0, 0], Boxed -> False, Axes -> None, Lighting -> False]

This animated sun effect was generated by translating the noise over time. Here is some sample POV-Ray code:
// runtime: 0 second
camera{location <0,0,2.25> look_at <0,0,0>}
sphere {<0,0,0>,1 pigment{bumps scale 0.1} finish{ambient 1}}

Julia Set Fractal - Mathematica 4.2: 4/15/04, POV-Ray 3.6.1: 7/4/06

z_n+1 = z_n²+z_c, z_c = -0.63-0.407i
Here is some Mathematica code for a typical Julia Set fractal. Julia Sets were discovered by Gaston Maurice Julia.
(* runtime: 17 seconds *)
Julia[zo_] := Module[{z = zo, i = 0}, While[i < 100 && Abs[z] < 2, z = z^2 + zc; i++]; i];
zc = -0.63 - 0.407 I;
DensityPlot[Sqrt[Julia[x + I y]], {x, -1.5, 1.5}, {y, -1.5, 1.5}, PlotPoints -> 275, Mesh -> False, Frame -> False, ColorFunction -> Hue]

POV-Ray has a built-in function for these fractals:
// runtime: 0.5 second
camera{orthographic location <0,0,-1.5> look_at 0 angle 90}
plane{z,0 pigment{julia <-0.63,-0.407>,200} finish{ambient 1}}

Cubic Julia Set Fractal - original version: Java, 6/9/01; animated version: Mathematica 4.2, 2/18/05

z_n+1 = z_n³+z_c, z_c = -0.5-0.05i
This was one of my older favorites. Here is some Mathematica code:
(* runtime: 10 seconds *)
Julia = Compile[{{z, _Complex}}, Length[FixedPointList[#^3 + zc &, z, 100, SameTest -> (Abs[#] > 2 &)]]];
zc = -0.5 - 0.05 I;
DensityPlot[Julia[x + I y], {x, -1.15, 1.15}, {y, -1.15, 1.15}, PlotPoints -> 275, Mesh -> False, Frame -> False, ColorFunction -> (RGBColor[Min[2#^2, 1], Max[2# - 1, 0], Min[2#, 1]] &)]

This is how you can make this fractal in POV-Ray:
// runtime: 0.5 second
camera{orthographic location <0,0,-1.15> look_at 0 angle 90}
plane{z,0 pigment{julia <-0.5,-0.05>,30 exponent 3 color_map{[0 rgb 0][1/3 rgb <0,0,1>][2/3 rgb <1,0,1>][1 rgb 1]}} finish{ambient 1}}

Inverse Julia Set Fractal - Mathematica 4.2, 6/26/04

Here is the same fractal using the inverse Julia set technique. This technique is good for showing the edges of the set, but some regions are faint because they attract much slower. Here is some Mathematica code:
(* runtime: 4 seconds *)
InvertList[zlist_] := Flatten[Map[(# - zc)^(1/3){1, -(-1)^(1/3) , (-1)^(2/3)} &, zlist], 1];
zc = -0.5 - 0.05 I; zlist = {0}; Do[zlist = InvertList[zlist], {10}];
ListPlot[{Re[#], Im[#]} & /@ zlist, PlotStyle -> PointSize[0.005], AspectRatio -> 1, Axes -> None]

This can be improved using pruning. The following Mathematica code was adapted from Mark McClure’s Mathematica code:
(* runtime: 4 seconds *)
InvertList[zlist_] := Flatten[Map[Floor[130(# - zc)^(1/3){1, -(-1)^(1/3), (-1)^(2/3)}]/130 &, zlist], 1];
zc = -0.5 - 0.05 I; zlist = Nest[InvertList, {0}, 8]; zlist2 = zlist;
While[zlist2 =!= {}, zlist2 = Complement[InvertList[zlist2],zlist]; zlist = Union[zlist2, zlist]];
ListPlot[{Re[#], Im[#]} & /@ zlist,PlotStyle -> PointSize[0.005], AspectRatio -> 1, Axes -> None]

Belousov-Zhabotinsky (BZ) Reaction - Mathematica 4.2, 5/10/07

This periodic cellular automaton simulates a chemical reaction called a Belousov-Zhabotinsky reaction. Here is some Mathematica code:
(* runtime: 27 seconds *)
n = 50; SeedRandom[0]; image = Table[Random[Integer, {0, 255}], {n}, {n}];
Do[image = Table[z = image[[i, j]]; zlist = Flatten[image[[Mod[i - {0, 1, 2},n] + 1, Mod[j - {0, 1, 2}, n] + 1]]]; If[z == 255, 0, a = 9 - Count[zlist, 0]; If[z ==0, b = Count[zlist, 255]; Floor[(a - b)/3] + Floor[b/3], Min[Floor[Plus @@ zlist/a] + 45, 255]]], {i, 1, n}, {j, 1, n}]; ListDensityPlot[image, Mesh -> False, Frame -> False], {200}];

Links
Mathematica notebook - by Takashi Yoshino
Five Cellular Automata program - by Hermetic Systems

Activator-Inhibitor, C++, Mathematica 4.2, 5/24/07

This periodic cellular automaton can be used to simulate patterns that resemble leopard spots. Here is some Mathematica code:
(* runtime: 7 minutes *)
n = 80; w = 0.35; SeedRandom[0]; image = Table[Random[Integer, {0, 1}], {n}, {n}];
Do[image = Table[ID = AD = 0; Do[If[image[[Mod[i + di - 1, n] + 1, Mod[j + dj - 1, n] + 1]] == 1, r = Sqrt[di^2 + dj^2]; If[3 <= r <= 6, ID++]; If[r <= 3, AD++]], {di, -6, 6}, {dj, -6, 6}]; a = AD - w ID; Which[a > 0, 1, a < 0, 0, True, image[[i, j]]], {i, 1, n}, {j, 1, n}]; ListDensityPlot[image, Mesh -> False, Frame -> False], {10}];

Links
Activator-Inhibitor - Wolfram Demonstrations Project
Heated metal alloy simulation - by Thomas Wanner

Apollonian Gasket - Mathematica 4.2, 7/20/05

Here are four circles inverted about each other. This is kind of similar to a Wada Basin but not exactly. Notice the small Poincar� hyperbolic tilings throughout the picture. Click here to download a Mathematica notebook for this animation. Click here to download some POV-Ray code. Here is some Mathematica code:
(* runtime: 0.3 second *)
chains = {N[Append[Table[{{2E^(I theta), Sqrt[3]}}, {theta, Pi/6, 9Pi/6, 2Pi/3}], {{0, 2 - Sqrt[3]}}]]};
Reflect[{z2_, r2_}, {z1_, r1_}] := Module[{a = r1^2/((z2 -z1)Conjugate[z2 - z1] - r2^2)}, {z1 + a (z2 - z1), a r2}];
Do[chains = Append[chains, Table[Map[Reflect[#, chains[[1, j, 1]]] &, Flatten[Delete[chains[[i]], j], 1]], {j, 1, 4}]], {i, 1, 6}];
Show[Graphics[Table[{GrayLevel[i/7], Map[Disk[{Re[#[[1]]], Im[#[[1]]]}, #[[2]]] &, chains[[i]], {2}]}, {i, 1, 7}], AspectRatio -> 1, PlotRange -> {{-1, 1}, {-1, 1}}]]

Links
Circle Inversion Java applet - by Pat Hooper
Apollonian Gasket - by Paul Bourke
Apollonian movie - by David Dumas

Barnsley’s Fern - Mathematica 4.2, 10/16/04

This is an example of a Directed Graph Iterated Function System (Digraph IFS). The image is generated by randomly selecting transformations which jump to different points on the fern. Here is some Mathematica code:
(* runtime: 7 seconds *)
n = 275; p = {0, 0}; image = Table[0, {n}, {n}];
Do[x = Random[Integer, 99]; p = Which[x < 3, {{0, 0}, {0, 0.16}}.p,x < 76, {{0.85, 0.04}, {-0.04, 0.85}}.p + {0, 1.6}, x <89, {{0.2, -0.26}, {0.23, 0.22}}.p + {0, 1.6}, True, {{-0.15, 0.28}, {0.26, 0.24}}.p + {0, 0.44}]; image[[Floor[n(p[[1]]/10 + 0.5)] + 1, Floor[n p[[2]]/10] + 1]]++, {100000}];
Show[Graphics[RasterArray[Map[RGBColor[0, 1 - Exp[-0.1#], 0] &, image, {2}]], AspectRatio -> 1]]

Links
Video Feedback - amazing video feedback fractals by Peter King
Mathematica code - by Mark McClure

Flame Fractal - Mathematica 4.2, 3/5/05

Flame Fractals are a popular, generalized variation of IFS fractals, invented by Scott Draves. This one is based on the Apophysis sample parameters. See also my POV-Ray code. Here is some Mathematica code:
(* runtime: 25 seconds *)
V1[{x_, y_}] := {Sin[x], Sin[y]};
F1[p_] := Module[{p2 = {{0.81, -0.33}, {0.08, 0.89}}.p + {0.24, -0.07}}, 0.88p2 + 0.12V1[p2]];
F2[p_] := Module[{p2 = {{0.3, 0.52}, {-0.56, 0.37}}.p + {-0.09, -0.42}}, 0.88p2 + 0.12V1[p2]];
F3[p_] := V1[{{-0.43, 0.38}, {-0.2, -0.44}}.p + {1.74, 1.21}];
F4[p_] := V1[{{-0.49, -0.27}, {0.28, -0.53}}.p + {0.51, 0.62}];
w1 = 0.65; w2 = 0.29; w3 = 0.03; w4 = 0.03;
n = 275; image = Table[{0, 0, 0}, {n}, {n}]; p = {0, 0};
Do[x = Random[]; k = Which[x < w1, 1, x < w1 + w2, 2, x < w1 + w2 + w3, 3, True, 4]; p = {F1, F2, F3, F4}[[k]][p]; {j,i} = Round[n{p[[1]] + 1.1, p[[2]] + 1.8}/3]; image[[i, j]] += List @@ ToColor[Hue[k/4.0], RGBColor], {100000}];
Show[Graphics[RasterArray[Map[RGBColor @@ Map[(1 - Exp[-0.2#]) &, #] &, image, {2}]], AspectRatio -> 1]]

Links
Cory Ench - some impressive flame fractals
Apophysis - free Flame Fractal software
Flame Fractal Java applet

Clifford Attractor - adapted from Paul Richards, Mathematica 4.2, 12/27/04

x_n+1 = sin(a y_n) + c cos(a x_n)
y_n+1 = sin(b x_n) + d cos(b y_n)
This type of strange attractor was invented by Clifford Pickover. If you look closely, you may find what appears to be distorted bifurcation fractals in the image. The right image is a superposition of many Clifford Attractors. See also my POV-Ray code. Here is some Mathematica code:
(* runtime: 3 minutes *)
n = 275; image = Table[{0, 0, 0}, {n}, {n}];
Interpolate[p_, x1_, x2_] := 2Cos[ArcCos[x1/2] + p(ArcCos[x2/2] - ArcCos[x1/2])];
Do[a = Interpolate[p,1.6, 1.3]; b = Interpolate[p, -0.6, 1.7]; c = Interpolate[p, -1.2, 0.5]; d = Interpolate[p, 1.6, 1.4]; x = y = 0; Do[{x, y} = {Sin[a y] + c Cos[a x], Sin[b x] + d Cos[b y]}; {i, j} = Round[n({y, x}/6 + 0.5)]; If[0 < i <= n && 0 < j <= n,image[[i, j]] += List @@ ToColor[Hue[p], RGBColor]], {1000}], {p, 0, 1, 0.001}];
Show[Graphics[RasterArray[Map[RGBColor @@ Map[1 - Exp[-0.02#] &, List @@ #] &, image, {2}]], AspectRatio -> 1]]

Hénon Map Escape Time - as seen on MathWorld, Mathematica 4.2, 6/16/04

x_n+1 = 1 - a x_n² + y_n, a = 0.2
y_n+1 = b x_n, b = 1.01
This beautiful fractal resembles a turbulent vortex. I saw this beautiful image on MathWorld but the image quality was poor due to aliasing. All I have done here is reproduce the image with antialiasing. The Hénon Map is defined by the above equations. Here is some Mathematica code:
(* runtime: 12 minutes *)
Henon[{x_, y_}] := {1 - 0.2 x^2 + y, 1.01 x};
DensityPlot[Length[FixedPointList[Henon[#] &, {x, y}, 1000, SameTest -> (#[[1]]^2 + #[[2]]^2 > 100 &)]], {x, 0, 4}, {y, -4, 0}, PlotPoints -> 275, AspectRatio -> 1, Mesh -> False, Frame -> False, ColorFunction -> (Hue[1 - #] &)]

Lorenz Attractor - Mathematica 4.2, 4/10-/06

The Lorenz Attractor is a famous chaotic strange attractor. The equations were originally developed to model atmospheric convection. Here is some Mathematica code:
(* runtime: 1 second *)
sigma = 3; rho = 26.5; beta = 1;
soln = {x[t], y[t], z[t]} /. NDSolve[{x'[t] == sigma (y[t] - x[t]), y'[t] == rho x[t] - x[t]z[t] - y[t], z'[t] == x[t] y[t] - beta z[t], x[0] == 0, y[0] == 1, z[0] == 1}, {x[t], y[t], z[t]}, {t, 0, 100}, MaxSteps -> 10000][[1]];
ParametricPlot3D[soln, {t, 0, 100}, PlotPoints -> 10000, Compiled -> False]

Here is some Mathematica code to numerically solve this using the 4th order Runge-Kutta method. See also my POV-Ray code:
(* runtime: 0.2 second *)
sigma = 3; rho = 26.5; beta = 1; dt = 0.01; p = {0, 1, 1}; v[{x_, y_, z_}] = {sigma(y - x), rho x - x z - y, x y - beta z};
Show[Graphics3D[Line[Table[k1 = dt v[p]; k2 = dt v[p + k1/2]; k3 = dt v[p +k2/2]; k4 = dt v[p + k3]; p += (k1 + 2 k2 + 2 k3 + k4)/6, {1000}]]]]

Links
Chaoscope - beautiful strange attractor software by Nicolas Desprez
Lorenz System Fattened Trajectory - interesting manifold by Michael Henderson
Lorenz and Modular Flows - nice renderings by Jos Leys
Mathematica code - by Mark McClure
Java applet - by Patrick Worfolk

Crackle Fractal - Mathematica 4.2, 11/17/04

Here is some Mathematica code for generating a periodic cellular texture based on the Voronoi diagram. This pattern resembles giraffe skin. See also my Delaunay Triangulation:
(* runtime: 3.8 minutes *)
SeedRandom[0]; nodes = Table[Random[], {100}, {2}];
Norm[x_] := x.x; Wrap[p_] := Map[If[# > 0.5, 1 - #, #] &, p];
DensityPlot[Module[{rlist = Sort[Map[Norm[Wrap[Abs[# - {x, y}]]] &, nodes]]}, rlist[[2]] - rlist[[1]]], {x, 0, 1}, {y, 0, 1}, PlotPoints -> 275, Mesh -> False, Frame -> False];

Mathematica’s ComputationalGeometry package can be used to generate Voronoi diagrams:
(* runtime: 1.3 seconds *)
<< DiscreteMath`ComputationalGeometry`; SeedRandom[0]; nodes = Table[Random[], {100}, {2}]; DiagramPlot[nodes, PlotRange -> {{0, 1}, {0, 1}}]

POV-Ray has a built-in function for the crackle fractal:
// runtime: 0.5 second
camera{orthographic location <0,0,-4> look_at 0 angle 90}
plane{z,0 pigment{crackle color_map{[0 rgb 0][0.5 rgb <0,1,0>][1 rgb 1]}} finish{ambient 1}}

Links
Cellular Portrait - Java applet by Paul Chew, see also Golan Levin's Cellular Portraits
Cellular Texture Tutorial - by Jim Scott

3D Mandelbrot Set - Mathematica 4.2, MathGL3d, 10/17/04

This Mandelbrot fractal was interpolated using Bill Rood’s formula for continuous escape time (cet) as described in the book “The Colours of Infinity”. Here is some Mathematica code:
cet = n + log₂ln(R) - log₂ln|z|
(* runtime: 3 seconds *)
R = 6;
image = ParametricPlot3D[Module[{z = 0.0, i = 0}, While[i < 100 && Abs[z] < R^2, z = z^2 + xc + I yc; i++]; cet = If[i != 100, i + (Log[Log[R]] - Log[Log[Abs[z]]])/Log[2], 0]; {xc, yc, 0.5Min[0.1cet, 1], {EdgeForm[], SurfaceColor[Hue[1 - 0.1cet]]}}], {xc, -2.0, 1.0}, {yc, -1.5, 1.5}, PlotPoints -> 64, Boxed -> False, Axes -> False, DisplayFunction -> Identity];
<< MathGL3d`OpenGLViewer`;
MVShow3D[image, MVNewScene -> True];

Quasifuchsian Limit Set - Mathematica 4.2, POV-Ray 3.6.1, 4/15/06

Here is some Mathematica code for a simple Quasifuchsian limit set:
(* runtime: 15 seconds *)
ta = 1.87 + 0.1I; tb = 1.87 - 0.1I; tab = (ta tb + Sqrt[ta^2tb^2 - 4(ta^2 + tb^2)])/2; z0 = (tab - 2)tb/(tb tab - 2ta + 2I tab);
b = {{tb - 2I, tb}, {tb, tb + 2I}}/2; B = Inverse[b]; a = {{tab, (tab - 2)/z0}, {(tab + 2)z0, tab}}.B; A = Inverse[a];
Affine[{z1_, z2_}] := 0.01 Round[(z1/z2)/0.01]; Children[{z_, n_}] := {Affine[{a, b, A, B}[[#]].{z, 1}], #} & /@ Delete[Range[4], {3,4, 1, 2}[[n]]];
ListPlot[{Re[#[[1]]], Im[#[[1]]]} & /@ Nest[Union[Flatten[Children /@ #, 1]] &, Table[{Affine[{a, b, A, B}[[i]].{0, 1}],i}, {i, 1, 4}], 12], AspectRatio -> Automatic, PlotRange -> All]

Link: 3D Maskit Slices - by Kentaro Ito

Mandelbrot Set Biomorph - Mathematica 4.2, 2/24/05

Biomorphs are cross-shaped orbit traps composed of Pickover Stalks invented by Clifford Pickover. Here is some Mathematica code:
(* runtime: 33 seconds *)
f[z_] := Min[Abs[{Re[z], Im[z]}]]/0.03;
OrbitTrap[zc_] := Module[{z = 0, i = 0}, While[i < 100 && Abs[z] < 100 && (i < 2 || f[z] > 1), z = z^2 + zc; i++]; Hue[Arg[z]/(2Pi), Min[1, 2 f[z]], Max[0, Min[1, 2(1 - f[z])]]]];
Show[Graphics[RasterArray[Table[OrbitTrap[xc + I yc], {yc, -1.5, 1.5, 3/274}, {xc, -2.0, 1.0, 3/274}]], AspectRatio -> 1]]

Links
Biomorphs - by Paul Bourke
Dragon on a Pedestal - by Damien Jones

Burning Ship Fractal - as seen on Paul Bourke’s web site, Mathematica 4.2, 8/5/04

x_n+1 = x_n² - y_n² - x_c, y_n+1 = 2 |x_n y_n| - y_c
This is really weird. I have no idea why it looks the way it does. It’s amazing what kind of pictures you can get with the right equations. Here is some Mathematica code:
(* runtime: 5 minutes *)
BurningShip[zc_] := Length[NestWhileList[{#[[1]]^2 - #[[2]]^2, 2Abs[#[[1]]#[[2]]]} - zc &, {0, 0}, #[[1]]^2 + #[[2]]^2 < 200 &, 1, 250]];
DensityPlot[-BurningShip[{xc, yc}], {xc, 1.71, 1.8}, {yc, -0.01, 0.08}, PlotPoints -> 275, Mesh -> False, Frame -> False]

Magnet Fractals - Mathematica 4.2: 6/26/04, POV-Ray 3.6.1: 7/4/06

These fractals were originally designed for predicting magnetic phase-transitions. Here is some Mathematica code:
Mandelbrot[zc_] := Length[FixedPointList[f[#, zc] &, 0, 100, SameTest -> (Abs[#] > 2 &)]];
Magnet[] := DensityPlot[Log[Mandelbrot[xc + I yc]], {xc, -1, 3}, {yc, -2, 2}, PlotPoints -> 275, Mesh -> False, Frame -> False, ColorFunction -> (If[# != 1, Hue[#], RGBColor[0, 0, 0]] &)];

(* Magnet 1: runtime: 5 minutes *)
f[z_, zc_] := (z^2 + zc - 1.0)^2 / (2 z + zc - 2)^2;
Magnet[]

(* Magnet 2: runtime: 15 minutes *)
f[z_, zc_] := (z^3 + 3 (zc - 1) z + (zc - 1) (zc - 2))^2 / (3 z^2 + 3 (zc - 2) z + (zc - 1) (zc - 2) + 1.0)^2;
Magnet[]

POV-Ray has a built-in function for these fractals:
//Magnet 1: runtime: 0.5 second
camera{orthographic location <1.25,0,-2> look_at <1.25,0,0> angle 90}
plane{z,0 pigment{magnet 1 mandel 50 interior 1,1 color_map{[0 rgb 0][1/6 rgb <0,0,1>][1/3 rgb 1]}} finish{ambient 1}}

//Magnet 2: runtime: 0 seconds
camera{orthographic location <1,0,-2> look_at <1,0,0> angle 90}
plane{z,0 pigment{magnet 2 mandel 50 color_map{[0 rgb 0][1/6 rgb <1,0,0>][1/3 rgb 1][1 rgb 1][1 rgb 0]}} finish{ambient 1}}

Fractal Crown - as described on Paul Bourke’s web site, equations by Roger Bagula - Mathematica 4.2, 8/13/04

Here is some Mathematica code for this fractal:
(* runtime: 10 seconds *)
n = 275; a = 5.0; b = Log[2]/Log[3]; image = Table[0, {n}, {n}];
Do[w = Sum[E^(I (-a)^k t)/a^(b k), {k, 1, 14}]; {i, j} = Floor[n({Re[w], Im[w]}/1.25 + 0.5)]; image[[i, j]] = Abs[w], {t, -Pi, Pi, 0.001}];
ListDensityPlot[image, Mesh -> False, Frame -> False]

This is how you can make this fractal in POV-Ray:
// runtime: 2 seconds
camera{location <0,0,1.25> look_at 0}
#macro Exp(Z) exp(Z.x)* #end
#declare a=5; #declare b=ln(2)/ln(3); #declare t1=-pi;
#while(t1*pow(-a,k)*t1)/pow(a,b*k); #declare k=k+1; #end sphere{w,0.003 pigment{rgb vlength(w)} finish{ambient 1}} #declare t1=t1+0.001; #end

Tree Fractal - adapted from from Renan Cabrera’s Mathematica code, Mathematica 4.2, 10/17/04

This is an example of a bracketed Lindenmayer System (L-System). You can also use this technique to make lightning bolts. Here is some Mathematica code:
(* runtime: 1 second *)
Polar[{x_, y_}, theta_, r_] := {x + r Cos[theta], y + r Sin[theta]};
branch[x_] := x; branch[tree_List] := branch /@ tree;
branch[Line[{p1_, p2_}]] := Module[{r = Sqrt[(p2[[1]] - p1[[1]])^2 + (p2[[2]] - p1[[2]])^2], theta = ArcTan[p2[[1]] - p1[[1]], p2[[2]] - p1[[2]]]}, {Thickness[0.05r], RGBColor[0.5r, 1 - 0.75r, 0], Line[{p2, Polar[p2, theta - Pi/6, 0.8r]}], Line[{p2, Polar[p2, theta + Pi/4, 0.7r]}]}];
Show[Graphics[{Thickness[0.07], RGBColor[0.5, 0.25, 0], NestList[branch, Line[{{0, 0}, {0, 1}}], 12]}], AspectRatio -> 1, Background -> RGBColor[0, 0, 0], PlotRange -> {{-3, 3}, {0, 6}}]

Links
3D Tree Fractals - POV-Ray source code by Oliver Knill
Tree Macros - more POV-Ray source code
Lace - interesting flower fractal Java applet by Chris Laurel

Bifurcation Diagram (Feigenbaum Fractal) - Mathematica 4.2, 10/17/02

This fractal is based on the logistic equation:
x_n+1 = a x_n (1 - x_n)
Here is some Mathematica code:
(* runtime: 1 minute *)
Logistic[x_] := a x (1 - x); plist = {};
Do[x = 0.25; Do[x = Logistic[x], {100}]; Do[x = Logistic[x]; plist = Append[plist, {a, x}], {100}], {a, 2.8, 4, 10^-2}];
ListPlot[plist]

Links
Feigenbaum Function - a fractal based on the bifurcation diagram
Bifurcation Diagram Java Applet
Similarities to the Mandelbrot Set

Koch’s Snowflake - AutoCAD, AutoLisp, 10/25/00

This is another type of L-System.

Links
The self-similar rippling of leaf edges and torn plastic sheets - hyperbolic surfaces resembling Koch’s snowflake
L-Systems - Mathematica code by unknown author, hosted by Jan Poland

Web - Java applet, 5/22/01

This was the first “complex” fractal I ever made. I accidentally made this fractal while trying to generate the Mandelbrot Set. Unfortunately, I don’t remember how to make it. Please let me know if you figure it out! I think the original code went something like this:
(* runtime: 38 seconds *)
Web[xc_, yc_] := Module[{x = 0, y = 0, i = 0}, While[i < 100 && (x^2 + y^2) < 4, x = x^2 - y^2 + xc; y = 2x y + yc; i++]; i];
DensityPlot[Log[Web[xc, yc]], {xc, -2, 1}, {yc, -1.5, 1.5}, PlotPoints -> 275, Mesh -> False, Frame -> False, ColorFunction -> (If[# == 1, Hue[0, 0, 0], Hue[0,Min[1, 2(1 - #)], Min[1, 2#]]] &)]

Want to see more? There were too many fractals on this page so I moved some. Click here to see more fractals along with Mathematica code:

More Fractals
Fractals Links
The Art of Mathematics - fractal slideshow narrated by Lasse Rempe
Fractal Art Contests - managed by Damien Jones
Jock Cooper - another good fractal artist
String Fractals - by Jos Leys
Summary of Fractal Types - a long list by Noel Giffin
Fractal Java applets
Mandelbrot & Julia Set Java Applet - by Mark McClure