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.
The Mandelbrot set was discovered in 1980 by Benoît Mandelbrot and is the most famous of all fractals. It is defined by iterating the function f(z) = z2 + c. For example, the third level Mandelbrot polynomial is given by F3(z) = f(f(f(z))). The Mandelbrot set is usually visualized using the Escape Time Algorithm (ETA) but another unique way to visualize this fractal is by its orbits, which are shaped like a bunch of circles and cardioids. The nth cycle orbits are the regions of c values where |Fn'(zo)| < 1 and zo is any root of Fn(zo) = zo (an attractive fixed point). The following Mathematica code was adapted from Mark McClure'sMathematica notebook:
(* runtime: 2.7 seconds *)
F[z0_] := Nest[#^2 + c &, z0, n]; ToPt[z_] := {Re[z], Im[z]}; used = {};
Show[Graphics[Table[roots = Select[c /. NSolve[F[0] == 0, c], ! MemberQ[Chop[used - #], 0] &]; used = Join[used,roots]; {Hue[n/5], Map[Polygon[Table[ToPt[c /. FindRoot[{F[z] == z, D[F[z], z] == Exp[I theta]}, {c, #}, {z, #}]], {theta, 0, 2Pi, Pi/18}]] &, roots]}, {n, 1, 5}], AspectRatio -> Automatic]]
If only we could find some way to spread these pearls around in 3D, then perhaps we could create the elusive and long sought after "true 3D Mandelbrot set". A simple rotation around the x-axis looks kind of pretty, but it's certainly no true 3D Mandelbrot set. Instead, we need to find a way to rotate the pearls around their parent pearls with multiple branching and variable sizes. See also my 3D Mandelbulb orbit.
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]
zn+1 = zn2+c
This Mandelbrot set zoom animation was created using the popular Escape Time Algorithm (ETA). The animation zooms through 15 orders of magnitude. Double precision numbers are inadequate for such high resolutions, so I had to use “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[c_] := Module[{z = 0, i = 0}, While[i < 100 && Abs[z] < 2, z = z^2 + c; 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]] &)]
This beautiful variation of the Mandelbrot set kind of resembles a nebula or a hologram. The Nebulabrot is a variation of the Buddhabrot technique. The coloring method that I used here was adapted from Paul Bourke’s web site. See also my 3D Nebulabrots. Here is some Mathematica code:
(* runtime: 10 minutes *)
Mandelbrot[c_] := Length[NestWhileList[#^2 + c &, 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]]
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]
Chop[f[z]]
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]
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]
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[{c = x + I y, z = 0.0}, Do[z = z^2 + c, {25}]; RGBColor @@ (image[[Floor[n Mod[Re[Log[Log[z]]]m/n, 1]] + 1, Floor[m Mod[2Arg[c], 1]] + 1]]/255.0)], {y, -1.5, 1.5, 3/275}, {x, -2, 1, 3.0/275}]], AspectRatio -> 1]]
zn+1 = zn2+c, c = -0.63-0.407i
Julia sets were discovered by Gaston Maurice Julia. Here is some Mathematica code:
(* runtime: 17 seconds *)
Julia[zo_] := Module[{z = zo, i = 0}, While[i < 100 && Abs[z] < 2, z = z^2 + c; i++]; i];
c = -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}}
zn+1 = zn3+c, c = -0.5-0.05i
This was one of my older favorites. Here is some Mathematica code demonstrating a simple technique called the Escape Time Algorithm (ETA):
(* runtime: 10 seconds *)
Julia = Compile[{{z, _Complex}}, Length[FixedPointList[#^3 + c &, z, 100, SameTest -> (Abs[#] > 2 &)]]];
c = -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}}
Here is another technique called the Cauchy Convergence Algorithm (CCA). This method uses orbit stability to determine whether or not a given point is inside the Julia set:
(* runtime: 9.4 minutes *)
f[z_] := z^3 + c; c = -0.5 - 0.05 I;
d[z_, w_] := 2ArcTan[Abs[(z Conjugate[w] - Abs[w]^2)/(z Conjugate[w] + 1)]/Abs[w]];
image = Table[z = x + I y; w = z + 1.0*^-10; d[z, w] + Sum[z = f[z]; w = f[w]; d[z, w], {18}], {y, -1.5, 1.5, 3.0/275}, {x, -1.5, 1.5, 3.0/275}];
ListDensityPlot[Log[image], Mesh -> False, Frame -> False]
The following Mathematica code demonstrates how to create a Julia set using the Inverse Iteration Method (IIM). This technique is good for showing the edges of the set, but some regions are faint because they attract much slower:
(* runtime: 4 seconds *)
c = -0.5 - 0.05 I; zlist = {0}; Do[zlist = Flatten[Map[(# - c)^(1/3){1, -(-1)^(1/3) , (-1)^(2/3)} &, zlist], 1], {10}];
ListPlot[{Re[#], Im[#]} & /@ zlist, PlotStyle -> PointSize[0.005], AspectRatio -> Automatic]
There are a variety of Modified Inverse Iteration Methods (MIIM) which can be used to speed up the calculation by pruning dense branches of the tree. One such method is to prune branches when the map becomes contractive (the cumulative derivative becomes large). The following Mathematica code was adapted from Peter Liepa'ssource code:
(* runtime: 3 seconds *)
power = 3; c = -0.5 - 0.05 I; zlist = {}; imax = 100; z = Table[0, {imax}]; dz = Table[1, {imax}]; roots = Table[1, {imax}]; i = 2;
While[i > 1, z[[i]] = {1, -(-1)^(1/3), (-1)^(2/3)}[[roots[[i]]]](z[[i - 1]] - c)^(1/power); dz[[i]] = power Abs[z[[i]]]^(power - 1)dz[[i - 1]]; zlist =Append[zlist, z[[i]]]; If[i < imax && dz[[i]] < 250.0, i++; roots[[i]] = 1, While[i > 1 && roots[[i]] == power, roots[[i]] = 1; i--]; roots[[i]]++]];
ListPlot[{Re[#],Im[#]} & /@ zlist, PlotStyle -> PointSize[0.005], AspectRatio -> Automatic]
Here is another approach for creating MIIM Julia sets, based on pixel hit counts:
(* runtime: * minutes *)
Here is another approach for creating MIIM Julia sets, adapted from Mark McClure’sMathematica code:
(* runtime: 4 seconds *)
InvertList[zlist_] := Flatten[Map[Floor[130(# - c)^(1/3){1, -(-1)^(1/3), (-1)^(2/3)}]/130 &, zlist], 1];
c = -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 -> Automatic]
The Mandelbrot set is defined by iterating the function f(z) = z2 + c. For example, the third level Mandelbrot polynomial is given by F3(z) = f(f(f(z))). Now, suppose we arbitrarily choose a "seed" value of z = 0 and solve for the roots of c such that Fn(0) = 0. Here is a plot of the roots of the 12th level Mandelbrot polynomial. Altogether, there are 2048 roots. Of course, this is still a far cry from the inverse Julia set pruning methods which can probe hundreds or even thousands of levels deep. Unfortunately, I can't think of any similar way to prune Mandelbrot set roots. Here is some Mathematica code:
(* runtime: 30 seconds *)
imax = 8; z = c /. Solve[Nest[#^2 + c &, 0.0, imax] == 0, c];
ListPlot[Map[{Re[#], Im[#]} &, z], AspectRatio -> 1]
You can also find these roots using the Durand-Kerner method:
(* runtime: 38 seconds *)
imax = 7; order = 2^(imax - 1); F[c_] := Nest[#^2 + c &, 0.0, imax]; z = Table[(0.4 + 0.9I)^i, {i, 1, order}]; z0 = 0.0z;
While[Max[Abs[z - z0]] > 1.0*^-7, z0 = z; z -= F[z]/Table[Product[If[i != j, z[[i]] - z[[j]], 1], {j, 1, order}], {i, 1, order}]];
ListPlot[Map[{Re[#], Im[#]} &, z], AspectRatio -> 1]
The Glynn Set is a special Julia Set, named after Earl Glynn. It kind of looks like a tree. See also my 3D Glynn Julia set. Here is some Mathematica code:
(* runtime: 48 seconds *)
Julia[z0_] := Module[{z = z0, i = 0}, While[i < 100 && Abs[z] < 2, z = z^1.5 + c; i++]; i];
c = -0.2;
DensityPlot[Julia[-x + I y], {y, -0.2, 0.2}, {x, 0.35, 0.75}, PlotPoints -> 275, Mesh -> False, Frame -> False]
Here is some code to plot this using the Modified Inverse Iteration Method (MIIM). Note that special care must be taken to verify each root's validity:
(* runtime 7 seconds *)
Pow[z_, n_, k_] := Module[{theta = Arg[z]}, theta = n(theta + 2Pi (k - Floor[(theta/Pi +Abs[1/n])/2])); If[Abs[theta] > Pi, Null, Abs[z]^n Exp[I theta]]];
power = 1.5; nroot = 3; c = -0.2; zlist = {}; dzmax = 25.0; imax =1000; z = Table[-0.61, {imax}]; dz = Table[1, {imax}]; roots = Table[1, {imax}]; i = 2;
While[i > 1, z[[i]] = Pow[z[[i - 1]] - c, 1/power, roots[[i]] - 1]; If[z[[i]] === Null, prune = True, dz[[i]] = Abs[power]Abs[z[[i]]]^(power - 1)dz[[i - 1]]; zlist = Append[zlist, z[[i]]]; prune = (i == imax || dz[[i]] > dzmax)]; If[prune, While[i > 1 && roots[[i]] == nroot, roots[[i]] = 1; i--]; roots[[i]]++, i++; roots[[i]] = 1]];
ListPlot[{Re[#], Im[#]} & /@ zlist, PlotStyle -> PointSize[0.005], AspectRatio ->Automatic]
These ones are incorrect, but I still think they look interesting.
zn+1 = (|xn| + i |yn|)2 - c = xn2 - yn2 + 2 i |xnyn| - c
This weird picture is a modified form of the Mandelbrot set. It’s amazing what kind of pictures you can get with the right equations. Here is some Mathematica code:
(* runtime: 5 minutes *)
BurningShip[c_] := Module[{x = 0, y = 0, i = 0}, While[x^2 - y^2 < 200 && i < 250, {x, y} = {x^2 - y^2, 2Abs[x y]} - c; i++]; i];
DensityPlot[-BurningShip[{xc, yc}], {xc, 1.71, 1.8}, {yc, -0.01, 0.08}, PlotPoints -> 275, Mesh -> False, Frame -> False]
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]]
xn+1 = sin(a yn) + c cos(a xn) yn+1 = sin(b xn) + d cos(b yn)
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]]
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] &)]
Here is some Mathematica code to make it look like granite:
(* runtime: 2 minutes *)
ListDensityPlot[Table[f = 2^i; Sum[Abs[0.5 - Noise[f x, f y]]/f, {i, 0, 5}], {x, 0.0, 1.0, 1.0/n}, {y, 0.0, 1.0, 1.0/n}], Mesh -> False, Frame -> False]
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}}
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}}]]
xn+1 = 1 - a xn2 + yn, a = 0.2 yn+1 = b xn, 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 - #] &)]
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}]]]]
This Mandelbrot set height field 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 + log2ln(R) - log2ln|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];
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