data >> "data.m"

The Julia set for the Hagihara-Hawkins function

Our objective is to generate the Julia set of the function

R(z)=(((z^3)-z)/(1+4 z-(z^2))).

This function is tricky because it has no attractive orbits. It does have fixed parabolic points at 0 and ∞, thus we can generate a Julia set using the Leau-Fatou/Big step technique.

Preliminaries

R(z) has a parabolic fixed point at zero:

R[z_]=(((z^3)-z)/(1+4*z-(z^2)));
{R[0],R'[0]}

Series[R[R[z]],{z,0,8}]

and a parabolic fixed point at ∞:

RInf[z_] = Simplify[1/R[1/z]];
{RInf[0],RInf'[0]}

Since both of these have multiplier -1, we will deal largely with the function (R^2)(z)=R(R(z)). Here are those function together with compiled versions:

R2[z_] = Simplify[R[R[z]]];
R2Inf[z_] = Simplify[1/R2[1/z]];
cR2 = With[{r2=R2[z]},Compile[{{z,_Complex}},r2]];
cR2Inf = With[{r2inf=R2Inf[z]},Compile[{{z,_Complex}},r2inf]];

Setting up series for big steps

Here we define the functions F[k][z], which approximate (R^(2^(k+1)))(z)

Remove[F];
F[0][z_] = N[Normal[Series[R2[z],{z,0,30}]]];
Do[
F[k][z_] = Chop[N[Normal[Series[F[k-1][F[k-1][z]],{z,0,30}]]],10^-8],
{k,1,10}]

And we do a similar thing near ∞.

Remove[FInf];
FInf[0][z_] = N[Normal[Series[R2Inf[z],{z,0,30}]]];
Do[
FInf[k][z_] = Chop[N[Normal[Series[FInf[k-1][FInf[k-1][z]],{z,0,30}]]],10^-8],
{k,1,10}]

It’s generally a good idea to check that these types of things work using code like the following:

Table[
{k,
Last[F[k][z]]/. z->1.0/2^(k-1),
t1 = Timing[F[k][1.0/2^(k-1)]],
t2= Timing[Nest[R2,1.0/2^(k-1),(2^k)]],
t2[[1]]/t1[[1]]},
{k,1,10}]

If we examine the output at k=5, we see that F[5][z] looks to be a good approximation of (R^6)(z) for z=1/(2^4) (and, hopefully, for |z|<1/(2^4)). The last value indicates how many times vaster using this should be. Things are even better for larger values of k. Here’s similar code near ∞.

Table[
{k,
Last[FInf[k][z]]/. z->1.0/2^(k),
t1 = Timing[FInf[k][1.0/2^(k)]],
t2= Timing[Nest[R2Inf,1.0/2^(k),(2^k)]],
t2[[1]]/t1[[1]]},
{k,1,10}]

Setting up the compiled iteration functions

The compiled function FF uses the magnitude of z to decide which function or approximation to use and applies it. The function iterate performs the iterateion. It’s compiled and listable to run in parallel.

FF = With[{
f5=Function[z,Evaluate[F[5][z]]],
f6=Function[z,Evaluate[F[6][z]]],
f7=Function[z,Evaluate[F[7][z]]],
f8=Function[z,Evaluate[F[8][z]]],
f9=Function[z,Evaluate[F[9][z]]],
f10 = Function[z,Evaluate[F[10][z]]],
fInf5=Function[z,Evaluate[FInf[5][z]]],
fInf6=Function[z,Evaluate[FInf[6][z]]],
fInf7=Function[z,Evaluate[FInf[7][z]]],
fInf8=Function[z,Evaluate[FInf[8][z]]],
fInf9=Function[z,Evaluate[FInf[9][z]]],
fInf10 = Function[z,Evaluate[FInf[10][z]]]
},
Compile[{{z,_Complex}},
Which[
1/(2^4)<Abs[z]<2^7,cR2[z], 
Abs[z]≤1/(2^9),f10[z],
Abs[z]≤1/(2^8),f9[z],
Abs[z]≤1/(2^7),f8[z],
Abs[z]≤1/(2^6),f7[z],
Abs[z]≤1/(2^5),f6[z],
Abs[z]≤1/(2^4),f5[z],
Abs[z]≥(2^9),1/fInf10[1/z],
Abs[z]≥(2^8),1/fInf9[1/z],
Abs[z]≥(2^7),1/fInf8[1/z],
Abs[z]≥(2^6),1/fInf7[1/z],
True,0],
RuntimeOptions -> "Speed", CompilationTarget -> "C"]];
iterate = With[{FF=FF},
Compile[{{z0,_Complex}},
Module[{z,i},
z=z0;
i = 0;
While[1/(2^9)<Abs[z]≤(2^9) && i++<10000,
z = FF[z]];
z],
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed", CompilationTarget -> "C"]];

Computing the iteration data

This is the most computationally intensive part. The fastest approach, which takes a little more than 20 minutes on my setup, is to run in parallel as follows:

800*1300

t=AbsoluteTime[];
data = iterate[Table[x+I*y + 0.0I,{y,8,0,-0.01},{x,-8,5,0.01}]];//AbsoluteTiming
AbsoluteTime[]-t

Unfortunately, that makes it hard to monitor the progress. Here’s a sequential approach that takes longer but allows you to monitor the progress.

data =Monitor[
 Table[ymonitor=y;iterate[x+I*y + 0.0I],{y,8,0,-0.1},{x,-8,5,0.1}],
ymonitor];//AbsoluteTiming

data =Monitor[
 Table[ymonitor=y;iterate[x+I*y + 0.0I],{y,1.2,0,-0.002},{x,-1.2,1.2,0.002}],
ymonitor];//AbsoluteTiming

5146/60.

Classifying

Now, we’ll want to classify the results according the Leau-Fatou Flower theorem. We really don’t even need the full blown version since, looking at the power series representations, it’s pretty easy to see (R^2) is weakly attracting along the real and imaginary axes and weakly repelling on axes rotated (45^°). (Subscript[R, ∞]^2), by contrast is weakly attracting along the real axis and weakly repelling along the imaginary.

Series[R2[z],{z,0,6}]

Series[R2Inf[z],{z,0,6}]

Taking this all into account, we might write the following classify function:

classify = Compile[{{z,_Complex}},
Module[{result,z2},
result = 0;
z2 = z*Exp[Pi*I/4];
Which[
Abs[z]<1/8,
Which[
-Pi<Arg[z2]<-Pi/2,result=1,
-Pi/2<Arg[z2]<0,result=2,
0<Arg[z2]<Pi/2,result=3,
Pi/2<Arg[z2]<Pi,result=4,
True,result=0],
Abs[z]>8,
Which[
Re[z]>0,result=5,
Re[z]<0,result=6,
True,result=0],
True,result=0];
result]];

Boundary scanning

Finally, we’ll want to boundary scan and plot the results.

kernel = {
{1,1,1},
{1,-8,1},
{1,1,1}
};
top = Map[classify,data,{2}];
bot = Reverse[Rest[top]] /. {2->4,4->2};
classified = Join[top,bot];
convolvedData = ListConvolve[kernel,classified];

pic = ArrayPlot[Sign[Abs[convolvedData]], Frame -> False]

pic = ArrayPlot[Sign[Abs[convolvedData]], Frame -> False, Axes->True,
DataRange->{{-2,2},{-2,2}}, Epilog->{Red,PointSize[Large],Point[{0,0}]}]

Rotate[%,Pi/2]

NotebookSave[]

[NotebookDirectory[]]

data = Get[FileNameJoin[{NotebookDirectory[],"data.m"}]];

data[[1,1]]

both = Table[
If[convolvedData[[i,j]]≠0,0,
classified[[i+1,j+1]]],
{i,1,Length[convolvedData]},{j,1,Length[First[convolvedData]]}];

ArrayPlot[both, Frame->False, PlotRangePadding->None,
DataRange -> {{-8,5},{-8,8}},Epilog->{PointSize[Large],Red,Point[{0,0}]},
ColorRules->{1->Darker@Green,2->Darker[Red],3->Lighter[Green],4->Lighter@Yellow,5->Darker@Blue,6->Lighter[Blue],0->Black}]

ArrayPlot[both, Frame->False, PlotRangePadding->None,
DataRange -> {{-8,5},{-8,8}},
ColorRules->{1->ColorData[61,3],2->ColorData[61,1],3->ColorData[61,6],4->ColorData[61,5],5->ColorData[61,2],6->ColorData[61,4],0->Black}]

Export["hhShaded.png", Out[62], ImageSize->{Automatic,1600}]

Dimensions[both]

Projecting onto the Riemann sphere

To place this on the Riemann sphere, we recompute the data using points that are more amenable for projection onto the sphere. We’ll also shade the basins a bit to help the result look a bit better.

step=Pi/600.;
pts=Table[Cot[phi/2] Exp[I*theta],
{phi,step,Pi-step,step},{theta,0,Pi-step,step}];

sData = iterate[pts];//AbsoluteTiming

right = Map[classify,sData,{2}];

classified = Map[Join[Reverse[Rest[#]]/. {2->4,4->2},#]&,right];

kernel = {
{1,1,1},
{1,-8,1},
{1,1,1}
};
convolvedData = ListConvolve[kernel,classified];

both = Table[
If[convolvedData[[i,j]]≠0,0,
classified[[i+1,j+1]]],
{i,1,Length[convolvedData]},{j,1,Length[First[convolvedData]]}];

texture = Image[ArrayPlot[both, Frame->False, PlotRangePadding->None,
ColorRules->{1->Red,2->Darker[Red],3->Red,4->Darker[Red],5->Gray,6->LightGray,0->Black}]];

texture = Image[ArrayPlot[both, Frame->False, PlotRangePadding->None,
ColorRules->{1->ColorData[61,3],2->ColorData[61,1],3->ColorData[61,6],4->ColorData[61,5],5->ColorData[61,2],6->ColorData[61,4],0->Black}]]

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[texture]],Lighting->"Neutral",Axes->False]

Export["texture.png",texture, ImageSize->1200]