How to draw Fractal images of iteration functions on the Riemann sphere?

Prof. McClure, in the work "M. McClure, Newton's method for complex polynomials. A preprint version of a “Mathematical graphics” column from Mathematica in Education and Research, pp. 1–15 (2006)", discusses how Mathematica can be applied to iteration functions for obtaining the basins of attraction (or their fractal images). Below, I provide his code for the fractal image of the polynomial $p(z)=z^3-1$:

p[z_] := z^3 - 1;
theRoots = z /. NSolve[p[z] == 0, z]
cp = Compile[{{z, _Complex}}, Evaluate[p[z]]];
n = Compile[{{z, _Complex}}, Evaluate[Simplify[z - p[z]/p'[z]]]];
bail = 150;
orbitData = Table[
NestWhileList[n, x + I y, Abs[cp[#]] > 0.01 &, 1, bail],
{y, -1, 1, 0.01}, {x, -1, 1, 0.01}
];
numRoots = Length[Union[theRoots]];
sameRootFunc = Compile[{{z, _Complex}}, Evaluate[Abs[3 p[z]/p'[z]]]];
whichRoot[orbit_] :=
Module[{i, z},
z = Last[orbit]; i = 1;
Scan[If[Abs[z - #] < sameRootFunc[z], Return[i], i++] &, theRoots];
If[i <= numRoots, {i, Length[orbit]}, None]
];
rootData = Map[whichRoot, orbitData, {2}];
colorList = {{cc, 0, 0}, {cc, cc, 0}, {0, 0, cc}};
cols = rootData /. {
{k_Integer, l_Integer} :> (colorList[[k]] /. cc -> (1 - l/(bail + 1))^8),
None -> {0, 0, 0}
};
Graphics[{Raster[cols]}]


My main question is here. He nicely obtained the fractal images on the complex plane, while it would be an interesting challenge to obtain these images on the Riemann sphere, e.g.

It seems the complex plane in this case has been replaced by a sphere, but how? I will be thankful if someone could revise the code given above for obtaining such beautiful fractal images on the Riemann sphere. Any tips and tricks will be fully appreciated as well.

-
Would you care to share where your example spherical projection came from? –  Mark McClure Dec 24 '12 at 13:04

As the other answers have shown, it's fairly easy to map an image onto a parametrized surface using textures. It can be a bit tricky, though, getting the image to mesh well with the transformation. J.M. hit on the crucial issue, namely that we compute the image using points that map to the sphere with minimal distortion. This answer is largely an expansion on his, although there are some differences and other ideas as well.

First, the article that Fazlollah refers to is some years old now, and the code can be improved in light of the many changes since V5, so let's start by showing how to generate regular Newton iteration images for general polynomials. Given a polynomial function $f(z)$, the following code computes the corresponding Newton's method iteration function $n$. It then defines the command limitInfo that iterates $n$ up to $50$ times from a starting point $z_0$ terminating when $|f(z)|$ is small and returning the last iterate and the number of iterates required for $|f(z)|$ to get small. It's compiled, listable and set to run in parallel, so it should be pretty fast. In addition to the function f, there are two numeric parameters to set, bail and r.

f = Function[z, z^3 - 1]; (* A very standard example *)
n = Function[z, Evaluate[Simplify[z - f[z]/f'[z]]]];
bail = 50;
(* bail is the number of iterates before bailing. *)
(* Doesn't have to be particularly large, *)
(* if there are only simple roots. *)
r = 0.01;
(* We assume that if |z-z0|<r, then we've *)
(* converged to the root z0. *)
limitInfo = With[{bail = bail, r = r, f = f, n = n},
Compile[{{z0, _Complex}},
Module[{z, cnt},
cnt = 0; z = z0;
While[Abs[f[z]] > r && cnt < bail,
z = n[z];
cnt = cnt + 1
];
{z, cnt}],
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
]];


Since the function is listable and runs in parallel, we can simply apply it to a table of data on one fell swoop.

step = 4/801; (* The denominator is essentially the resolution. *)
limitData = limitInfo[
Table[x + I*y, {y, 2, -2, -step}, {x, -2, 2, step}]]; // AbsoluteTiming

(* Out: {1.492716, Null} *)


Each element is a pair that indicates the limiting behavior and how long it took to get there.

limitData[[1, 1]]

(* Out: {-0.499835 + 0.866044 I, 5. + 0. I} *)


I guess we need a function that takes something like that and turns it into a color.

roots = z /. NSolve[f[z] == 0, z];
preColors = List @@@ Table[ColorData[61, k], {k, 1, Length[roots]}];
preColors = Append[preColors, {0.0, 0.0, 0.0}];
color = With[{bail = bail, roots = roots, preColors = preColors},
Compile[{{z, _Complex}, {cnt, _Complex}},
Module[{arg, time, i},
arg = Arg[z];
time = Abs[cnt];
i = 1;
Scan[If[Abs[z - #] < 0.1, Return[i], i++] &, roots];
Abs[preColors[[i]]*(cnt/bail)^(0.2)]
(* The exponent 0.2 adjusts the brightness of the image. *)]]
];


Now, we apply that function and generate the image.

colors = Apply[color, limitData, {2}];
Image[colors, ImageSize -> 2/step]


To map onto a sphere nicely, we'll discard the rectangular grid of points that we used above in favor of a collection of points that looks something like the following (although, we'll want higher resolution, of course):

step = Pi/12;
pts = Table[Cot[phi/2] Exp[I*theta],
{phi, step, Pi - step, step}, {theta, -Pi, Pi, step}];
ListPlot[{Re[#], Im[#]} & /@ Flatten[pts, 1],
AspectRatio -> Automatic, PlotRange -> All,
Epilog -> {Red, Circle[]}]


The expression $\cot(\phi/2) e^{i\theta}$ is the stereographic projection of a point expressed in spherical coordinates $(1,\phi,\theta)$ onto the plane. As a result, the corresponding points on the sphere are nicely distributed. Note, for example, that the number of points inside and outside of the unit circle are the same.

Graphics3D[{{Opacity[0.8], Sphere[]},
Point[Flatten[Table[{Cos[theta] Sin[phi], Sin[theta] Sin[phi], Cos[phi]},
{phi, step, Pi - step, step}, {theta, -Pi, Pi, step}], 1]]}]


Now, we increase the resolution and use the same limitInfo and color functions as before.

step = Pi/500;
limitData = limitInfo[Table[Cot[phi/2] Exp[I*theta],
{phi, step, Pi - step, step}, {theta, -Pi, Pi, step}]];
colors = Apply[color, limitData, {2}];
rect = Image[colors, ImageSize -> 4/step]


The image looks a bit different, but it's perfect for use as a spherical texture.

ParametricPlot3D[{Cos[theta] Sin[phi], Sin[theta] Sin[phi], Cos[phi]} ,
{theta, -Pi, Pi}, {phi, 0, Pi}, Mesh -> None, PlotPoints -> 100,
Boxed -> False, PlotStyle -> Texture[Show[rect]],
Lighting -> "Neutral", Axes -> False]


We can incorporate all of this into a Module.

newtonSphere[fIn_, var_, resolution_, bail_: 50, r_: 0.01] := Module[
{f, n, limitInfo, color, colors, roots, preColors, step, limitData, rect},

f = Function[var, fIn];
n = Function[var, Evaluate[Simplify[var - f[var]/f'[var]]]];
limitInfo = With[{bailLoc = bail, rLoc = r, fLoc = f, nLoc = n},
Compile[{{z0, _Complex}},
Module[{z, cnt},
cnt = 0; z = z0;
While[Abs[fLoc[z]] > rLoc && cnt < bailLoc,
z = nLoc[z];
cnt = cnt + 1
];
{z, cnt}],
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
]];
roots = z /. NSolve[f[z] == 0, z];
preColors = List @@@ Table[ColorData[61, k], {k, 1, Length[roots]}];
preColors = Append[preColors, {0.0, 0.0, 0.0}];
color = With[{bailLoc = bail, rootsLoc = roots, preColorsLoc = preColors},
Compile[{{z, _Complex}, {cnt, _Complex}},
Module[{arg, time, i},
arg = Arg[z];
time = Abs[cnt];
i = 1;
Scan[If[Abs[z - #] < 0.1, Return[i], i++] &, rootsLoc];
preColorsLoc[[i]]*(cnt/bailLoc)^(0.2)
]]
];
step = Pi/resolution;
limitData = limitInfo[
Table[Cot[phi/2] Exp[I*theta], {phi, step, Pi - step,
step}, {theta, -Pi, Pi, step}]];
colors = Apply[color, limitData, {2}];
rect = Image[colors, ImageSize -> 4/step];
ParametricPlot3D[{Cos[theta] Sin[phi], Sin[theta] Sin[phi],
Cos[phi]} ,
{theta, -Pi, Pi}, {phi, 0, Pi}, Mesh -> None, PlotPoints -> 100,
Boxed -> False, PlotStyle -> Texture[Show[rect]],
Lighting -> "Neutral", Axes -> False]
];


Now, if I had to guess, I'd say that example image in the original post was generated by a small perturbation of $z^8-z^2$.

newtonSphere[(2 z/3)^8 - (2 z/3)^2 + 1/10, z, 500]


Here are a few more examples.

pic1 = newtonSphere[z^2 - 1, z, 401];
SeedRandom[1];
pic2 = newtonSphere[Sum[RandomInteger[{-3, 5}] z^k, {k, 0, 8}], z, 400];
pic3 = newtonSphere[z^10 - z^5 - 1, z, 400, 200];
pic4 = newtonSphere[z^5 - z - 0.99, z, 400];
GraphicsGrid[{
{pic1, pic2},
{pic3, pic4}
}]


In the top row, we see that the result for quadratic polynomials is typically rather boring while that for a random degree 8 polynomial can be quite cool. On the bottom right, we see a black region. The color function is setup to default to black when none of the roots are detected. This can certainly happen; in fact the Newton iteration function for this example has an attractive orbit of period 6 leading to the quadratic like Julia set seen in the image. Sometimes black can occur simply because we didn't iterate enough, which is why I used the optional fourth argument for the image in the bottom left.

-
Now imagine those hanging from a fractal christmas tree. (+1) –  Jens Dec 24 '12 at 17:40

Here is my modest attempt, based on the formulae for stereographic projection in this Wikipedia entry (where the north pole corresponds to the point at infinity) and using a technique similar to the one in this answer:

newtonRaphson = Compile[{{n, _Integer}, {c, _Complex}},
Arg[FixedPoint[(# - (#^n - 1)/(n #^(n - 1))) &, c, 30]]]

tex = Image[DensityPlot[
newtonRaphson[3, Cot[ϕ/2] Exp[I θ]], {θ, -π, π}, {ϕ, 0, π},
AspectRatio -> Automatic,
ColorFunction -> (Which[# < .3, Red, # > .7, Yellow, True, Blue] &),
Frame -> False, ImagePadding -> None, PlotPoints -> 400,
PlotRange -> All, PlotRangePadding -> None],
ImageResolution -> 256];

(* yes, I know that I could have used SphericalPlot3D[]... *)
ParametricPlot3D[{Sin[ϕ] Cos[θ], Sin[ϕ] Sin[θ], Cos[ϕ]}, {θ, -π, π}, {ϕ, 0, π},
Axes -> None, Boxed -> False, Lighting -> "Neutral", Mesh -> None,
PlotStyle -> Texture[tex], TextureCoordinateFunction -> ({#4, #5} &)]


-
Thanks for your response. There are two problems. 1. How can one zoom in on a particluar place on this sphere without lowering the quality to observe the fractal behaviour of the method? 2. The "spcae size" of the output image? In fact, how one can save as the output fractal image with low "disk size" without lowering the quality in EPS format? For example, for $n=8$, its size is more than 2MB! –  Fazlollah Soleymani Nov 22 '12 at 21:16
You'll have to play with PlotPoints and ImageResolution on your own, of course... –  Guess who it is. Nov 22 '12 at 23:13
Thanks, the problem is that when we reduce PlotPoints or ImageResolution, the quality gets down dramatically. I am looking for a fast way to obtain high quality pics, with small space size, just like the one given in the question. Rasterize@... is a good choice, but it disables the feature of rotating the 3D pic, and also gets the quality lower. –  Fazlollah Soleymani Nov 23 '12 at 8:20
Well, I'd say this is one of those "no such thing as a free lunch" things. If you want to be able to zoom, you certainly need high resolution, which will demand more of your computer's resources... if you want small file sizes, you'll have to sacrifice quality somewhat. Scylla and Charybdis, you know... –  Guess who it is. Nov 23 '12 at 9:49
img2 = ImageCrop[Image[Graphics[{Raster[cols]}, PlotRangePadding -> 0,
ImagePadding -> 0, ImageMargins -> 0]], {343, 343}];

SphericalPlot3D[1 , {u, 0, Pi}, {v, 0, 2 Pi}, Mesh -> None,
TextureCoordinateFunction -> ({#1, #2} &),
PlotStyle -> Directive[Specularity[White, 10], Texture[img2]],
Lighting -> "Neutral", Axes -> False, ImageSize -> 500]


where img2 is cropped version of the 2D image in OP's question.

Check Texture for more examples.

-
Thanks for your reply, but there are big white circles at the middle of each basin. They should not be here. Please rotate your image, and then you will see the incomplete fractal image. Can you solve this drawback? –  Fazlollah Soleymani Nov 22 '12 at 20:59
img2 is the cropped version of your 2D image: img2=ImageCrop[ Image[Graphics[{Raster[cols]}, PlotRangePadding -> 0, ImagePadding -> 0, ImageMargins -> 0]], {343, 343}]. –  kglr Nov 22 '12 at 21:50
Now that I think about it, one could have directly produced an image instead of going through the Raster[] route: Image[cols]. –  Guess who it is. Nov 24 '12 at 14:12

I've decided to write a simplification+extension of Mark's routine as a separate answer. In particular, I wanted a routine that yields Riemann sphere fractals not only for Newton-Raphson, but also its higher-order generalizations (e.g. Halley's method).

I decided to use Kalantari's "basic iteration" family for the purpose. An $n$-th order member of the family looks like this:

$$x_{k+1}=x_k-f(x_k)\frac{\mathcal D_{n-1}(x_k)}{\mathcal D_n(x_k)}$$

where

$$\mathcal D_0(x_k)=1,\qquad\mathcal D_n(x_k)=\begin{vmatrix}f^\prime(x_k)&\tfrac{f^{\prime\prime}(x_k)}{2!}&\cdots&\tfrac{f^{(n-2)}(x_k)}{(n-2)!}&\tfrac{f^{(n-1)}(x_k)}{(n-1)!}\\f(x_k)&f^\prime(x_k)&\ddots&\vdots&\tfrac{f^{(n-2)}(x_k)}{(n-2)!}\\&f(x_k)&\ddots&\ddots&\vdots\\&&\ddots&\ddots&\vdots\\&&&f(x_k)&f^\prime(x_k)\end{vmatrix}$$

As noted in that paper, the basic family generalizes the Newton-Raphson iteration; $n=1$ corresponds to Newton-Raphson, while $n=2$ gives Halley's method. (Relatedly, see also Kalantari's work on polynomiography.)

Here's a routine for $\mathcal D_n(x)$:

iterdet[f_, x_, 0] := 1;
iterdet[f_, x_, n_Integer?Positive] := Det[ToeplitzMatrix[PadRight[{D[f, x], f}, n],
Table[SeriesCoefficient[Function[x, f]@\[FormalX], {\[FormalX], x, k}], {k, n}]]]


Here is the routine for generating the Riemann sphere fractals:

Options[rootFractalSphere] = {ColorFunction -> Automatic, ImageResolution -> 400,
MaxIterations -> 50, Order -> 1, Tolerance -> 0.01};

rootFractalSphere[fIn_, var_, opts : OptionsPattern[]] /; PolynomialQ[fIn, var] :=
Module[{γ = 0.2, bail, cf, colList, f, h, itFun, ord, roots, tex, tol},

f = Function[var, fIn];
ord = OptionValue[Order];
itFun = Function[var, var - Simplify[f[var] iterdet[f[var], var, ord - 1]/
iterdet[f[var], var, ord]] // Evaluate];

roots = var /. NSolve[f[var], var];
cf = OptionValue[ColorFunction];
If[cf === Automatic, cf = ColorData[61]];
colList = Append[Table[List @@ ColorConvert[cf[k], RGBColor], {k, Length[roots]}],
{0., 0., 0.}];

bail = OptionValue[MaxIterations]; tol = OptionValue[Tolerance];
makeColor = Compile[{{z0, _Complex}},
Module[{cnt = 0, i = 1, z},
z = FixedPoint[(++cnt; itFun[#]) &, z0, bail,
SameTest -> (Abs[f[#2]] < tol &)];
Scan[If[Abs[z - #] < 10 tol, Return[i], i++] &, roots];
Abs[colList[[i]] (cnt/bail)^γ]],
CompilationOptions -> {"InlineExternalDefinitions" -> True},
RuntimeAttributes -> {Listable}, RuntimeOptions -> "Speed"];

h = π/OptionValue[ImageResolution];
tex = DeveloperToPackedArray[makeColor[
Table[Cot[φ/2] Exp[I θ], {φ, h, π - h, h}, {θ, -π, π, h}]]];

ParametricPlot3D[{Cos[θ] Sin[φ], Sin[θ] Sin[φ], Cos[φ]}, {θ, -π, π}, {φ, 0, π},
Axes -> False, Boxed -> False, Lighting -> "Neutral", Mesh -> None,
PlotPoints -> 75, PlotStyle -> Texture[tex],
Evaluate[Sequence @@ FilterRules[{opts}, Options[Graphics3D]]]]]


Other notes:

• The compiled functions limitInfo[] and color[] have been merged into the single function makeColor[]. This function was not localized on purpose to allow its use even after executing rootFractalSphere[].

• Texture[] can directly accept an array of RGB triplets, so there is no need to use Image[] if these triplets are being generated directly by makeColor[].

Now, for some examples. The first two are Newton-Raphson fractals:

rootFractalSphere[z^3 - 1, z]


rootFractalSphere[(2 z/3)^8 - (2 z/3)^2 + 1/10, z]


Here is a fractal generated by Halley's method:

rootFractalSphere[(2 z/3)^8 - (2 z/3)^2 + 1/10, z, Order -> 2]


Finally, a fractal from a third order iteration:

rootFractalSphere[z^10 - z^5 - 1, z, ColorFunction -> ColorData[54],
MaxIterations -> 200, Order -> 3]
`

-
Thanks. An excellent answer. Just one note. There might be some diverging points (black areas) in the fractal picture for some test problems. On the other hand, we used a color correspond to a point in a sphere or a rectangular domain. So, the domain of working (I mean the mesh of points) are finite. So is it possible to count the number of diverging points? I mean it would be nice to have the percentage of diverging points for each fractal picture. Is it possible to cunt the number of diverging points in your implementation? –  Fazlollah Soleymani Apr 6 '13 at 12:25