Grid[{{pic1},{pic2}}]

Computing parabolic Julia sets

In this notebook, we present some code for computing Julia sets that contain a parabolic fixed point. Specifically, we’ll look at functions of the form

f(z)=z+a (z^(n+1))+O((z^(n+2)))

for, then, the origin is a parabolic fixed point.

The Leau-Fatou flower theorem

The Leau-Fatou flower theorem then tells us the nature of the dynamics in a neighborhood of the origin. Specifically, consider the two equations

n a (v^n)=1
n a (v^n)=-1.

any vector v that is a solution of the first equation is a repulsion vector, while any vector v that is a solution of the second equation is an attraction vector. The repulsion vectors partition a small neighborhood of the origin into regions containing invariant pedals. We can classify the Fatou set according to which of these pedals the iteration lands.

A first example

One particularly simple example is f(z)=z+(z^5). It’s pretty easy to see by direct inspection that this function is repulsive in the directions ±1 and ±i; let’s see how this comes out of the Leau-Fatou flower theorem, though. Note that we have n=4 and a=1. Thus, it’s just a matter of solving.

Clear[v];
v/.Solve[4(v^4)==1,v]

This allows us to write the following function to classify the nature of an orbit after a few iterates. Here’s a function that does that.

simpleClassify = Compile[{{z,_Complex}},
Module[{arg},
If[Abs[z]>2,0, (* Orbit escaped *)
arg = Arg[z];
Which[
0<arg<Pi/2,1,(* landed in Quadrant 1 *)
Pi/2<arg<Pi,2,(* Quadrant 2 *)
-Pi<arg<-Pi/2,3,(* Quadrant 3 *)
-Pi/2<arg<0,4,(* Quadrant 4 *)
True,5]]]];

To apply it, we first compute a table of values after iteration.

itValue = Compile[{{z0,_Complex}},
Module[{z,cnt, bail},
z=z0;
cnt = 0;
If[Abs[z]<1/2,bail=100,bail=10];
While[Abs[z]≤2 && cnt++<bail,
z = z+(z^5)];
z]];
stepSize = 0.001;
data = Table[itValue[x+I*y],{y,1,-1,-stepSize},
{x,-1,1,stepSize}];//Timing

Here’s an image of the classification of points.

classified = Map[simpleClassify,data,{2}];
ArrayPlot[classified]

Now, we use a little boundary scanning to get the actual Julia set.

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

A more general classifier

The Leua-Fatou flower theorem allows us to do this in a much more general way.

classifyArg = Compile[
{{z,_Complex},{z0,_Complex},{v,_Complex},{n,_Integer}},
Module[{check,check2},
check=n(Arg[((z-z0)/(-v))]+Pi)/(2Pi);
check2=Ceiling[check];
If[check==check2,0,check2]],
CompilationTarget -> "C",
RuntimeOptions->"Speed"];

Eleven fold symmetry

Clear[f];
f[z_] = z+Exp[2Pi*I/5]*(z^12);
JuliaSetPlot[N[f[z]],z]

As expected, we see 11 fold symmetry.

Clear[a,v,n];
n = 11;
series = Normal[Series[f[z],{z,0,n+1}]]
a = Coefficient[series,z,n+1];
v = Last[v /. NSolve[a*(v^n)==1,v]]

f=Compile[{{z,_Complex}},
z+Exp[2Pi*I/5]*(z^12)];
it = With[{f=f},
Compile[{{z0,_Complex}},
Module[{z,i, bail},
z=z0;
i = 0;
If[Abs[z]<1/2,bail = 1000,bail=10];
While[Abs[z]≤2 && i++<bail,
z = f[z]];
z],CompilationTarget -> "C",
RuntimeOptions->"Speed",
RuntimeAttributes->{Listable}, Parallelization->True]];
stepSize = 0.001;
pts = Table[x+I*y,{y,1.2,-1.2,-stepSize},
{x,-1.2,1.2,stepSize}];
data = it[pts];//AbsoluteTiming

classify = With[{classifyArg = classifyArg,v=v,n=n},
Compile[{{z,_Complex}},
Which[
Abs[z]>2,0,
True,classifyArg[z,0.0+0.0I,v,n]],
CompilationTarget -> "C",
RuntimeOptions->"Speed"]];
classified = Map[classify,data,{2}];
convolvedData = ListConvolve[kernel,classified];
ArrayPlot[Sign[convolvedData],Frame->False,
AspectRatio -> Automatic]

A neutral fixed point combined with an attractive fixed point

Clear[f];
f[z_] = z+(0.8+0.8I)(z^4)+(z^5);
JuliaSetPlot[f[z],z]

This function has a parabolic fixed point at the origin but it also has an attractive fixed point.

z/.Solve[f[z]==z,z]

Chop[f'/@ %]

Clear[a,v,n];
n = 3;
series = Normal[Series[f[z],{z,0,n+1}]]
a = Coefficient[series,z,n+1];
v = Last[v /. NSolve[a*(v^n)==1,v]]

f=Compile[{{z,_Complex}},
z+(0.8+0.8I)(z^4)+(z^5)];
it = With[{f=f},
Compile[{{z0,_Complex}},
Module[{z,i},
z=z0;
i = 0;
While[Abs[z+(0.8+0.8I)]>0.05 && Abs[z]≤2 && i++<1000,
z = f[z]];
z],CompilationTarget -> "C",
RuntimeOptions->"Speed",
RuntimeAttributes->{Listable}, Parallelization->True]];
stepSize= 0.002;
pts = Table[x+I*y,{y,0.8,-1.2,-stepSize},
{x,-1.2,0.8,stepSize}];
data = it[pts];//AbsoluteTiming

classify = With[{classifyArg = classifyArg,v=v,n=n},
Compile[{{z,_Complex}},
Which[
Abs[z]>2,0,
Abs[z+(0.8+0.8I)]<0.1,4,
True,classifyArg[z,0.0+0.0I,v,n]],
CompilationTarget -> "C",
RuntimeOptions->"Speed"]];
classified = Map[classify,data,{2}];
convolvedData = ListConvolve[kernel,classified];
ArrayPlot[Sign[convolvedData],Frame->False,
AspectRatio -> Automatic]

A higher order mulitplier

Clear[f];
f[z_] = Exp[2Pi*I*3/7]z+(z^2);
JuliaSetPlot[N[f[z]],z]

Now, the multiplier is (e^(2π i 3/7)), rather than just 1. The Leau-Fatou flower theorem requires the multiplier to be 1 but, if we compute (f^7), we get a function with multiplier 1 and the same Julia set. We can check this via expansion.

Chop[N[Collect[Expand[Nest[f,z,7]],z]]][[1;;5]]

Generally, the more efficient way to do this is with a series approximation.

n=7;
series = Chop[N@Normal[Series[Nest[f,z,n],{z,0,n+1}]]]

Clear[a,v]
a = Coefficient[series,z,n+1];
v = Last[v /. NSolve[n*a*(v^n)==1,v]]

f=Compile[{{z,_Complex}},
Exp[2Pi*I*3/7]z+(z^2)];
it = With[{f=f},
Compile[{{z0,_Complex}},
Module[{z,i, bail},
z=z0;
i = 0;
While[Abs[z]≤2 && i++<7777, (* Need a large mulitple of 7 *)
z = f[z]];
z],CompilationTarget -> "C",
RuntimeOptions->"Speed",
RuntimeAttributes->{Listable}, Parallelization->True]];

pts = Table[x+I*y,{y,0.8,-1.2,-0.0025},
{x,-1.1,2.1,0.0025}];
data = it[pts];//AbsoluteTiming

classify = With[{classifyArg = classifyArg,v=v,n=n},
Compile[{{z,_Complex}},
Which[
Abs[z]>2,0,
True,classifyArg[z,0.0+0.0I,v,n]],
CompilationTarget -> "C",
RuntimeOptions->"Speed"]];
classified = Map[classify,data,{2}];

kernel = {
{1,1,1},
{1,-8,1},
{1,1,1}
};
convolvedData = ListConvolve[kernel,classified];
ArrayPlot[Sign[convolvedData],Frame->False,
AspectRatio -> Automatic]