Parabolic dynamics

The Leau-Fatou, big step method

A fixed point (Subscript[z, 0]) of the map f is called parabolic if f'((Subscript[z, 0])) is a root of unity. We can extend this to a point with period n by considering the map (f^n). Images involving parabolic dynamics can be difficult to generate since things move very slowly near the parabolic point. This notebook presents some theory and practical code to deal with the situation.

SetDirectory[FileNameJoin[{ParentDirectory[NotebookDirectory[]],"images"}]]

An beginning example

For our first example, we consider f(z)=z+(z^5). It’s pretty easy to see that the origin is a parabolic fixed point. Let’s try to plot it as easily as possible.

pic = JuliaSetPlot[z-32z^5,z,Frame->False,
PlotRange->{{-0.6,0.6},{-0.6,0.6}},PlotRangePadding->None]

Kinda cool but, ultimately, inaccurate. It’s not easy to improve it by simply increasing the iteation count; we’ll need to take a different approach.

It turn out that after a few iterates, most points have either diverged to infinity or landed in one of the four lobes attached to the origin. Thus, an alternative approach is to classify an initial seed based on which of those happens. Here’s some simple code to do this:

iterate = Compile[{{z0,_Complex}},
Module[{z,cnt, bail},
z=z0;
cnt = 0;
bail = 100;
While[Abs[z]≤2 && cnt++<bail,
z = z-32z^5];
z]];
stepSize = 0.002;
data = Table[iterate[x+I*y],{y,0.6,-0.6,-stepSize},
{x,-0.6,0.6,stepSize}];//AbsoluteTiming

Now, we’ll write a function to quickly classify the resulting landing points.

simpleClassify = Compile[{{z,_Complex}},
Module[{arg,z2},
If[Abs[z]>2,0, (* Orbit escaped *)
z2 = Exp[Pi*I/4]*z;
arg = Arg[z2];
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]]]];
classified = Map[simpleClassify,data,{2}];
ArrayPlot[classified, Frame -> False, PlotRangePadding->None]

Finally, to get the Julia set itself, we’ll employ a boundary scanning technique.

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

Export["boundaryScan32.png", %, ImageSize->{800,800}]

The Leau-Fatou flower

The Leau-Fatou flower theorem states that, if f has the Taylor expansion

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

then the complex number v is a repulsion vector for f at the origin if n a (v^n)=1 and v is an attraction vector if n a (v^n)=-1.

If we apply this to the z+(z^5) example, then a=1 and n=4 so we get

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

Here’s a fairly general classifier function we can use based on this theorem.

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

A harder set of examples

Let’s consider the family of functions f(z)=(e^(2π i/n))z+(z^2). While not directly applicable, the Leau-Fatou flower theorem does apply to the function (f^n). Let’s start with n=7.

n=7;
f[z_] = Exp[2Pi*I/n]z+z^2;
series = Chop[N[Series[Nest[f,z,n],{z,0,n+1}]],10^-8];
Clear[v,vv];
a = Coefficient[series,z,n+1];
vv = v/.Solve[n*a*v^n==1,v];
v = First[vv]

We’ll write a classify function that uses classifyArg near zero but also accounts for the possibility of divergence.

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]]]];

We’ll iterate 1000 times and then see what happens.

iterate =With[{n=n},
 Compile[{{z0,_Complex}},
Module[{z,cnt, bail},
z=z0;
cnt = 0;
bail = 1000;
While[Abs[z]≤2 && cnt++<bail,
z = Exp[2Pi*I/n]z+z^2];
z]]];

We’ll use some symmetry to cut computation time in half.

center=First[Select[z/.NSolve[f[z]==0,z],Im[#]<0&]]/2

Do it.

data = Monitor[
Table[iterate[x+I*y],{y,Im[center]+1.2,Im[center],-0.002},
{x,Re[center]-1.2,Re[center]+1.2,0.002}],y];//AbsoluteTiming

Use the symmetry.

data = Join[data,Reverse[Rest[Reverse/@data]]];

And make the image.

kernel = {
{1,1,1},
{1,-8,1},
{1,1,1}
};
classified = Map[classify,data,{2}];
convolvedData = ListConvolve[kernel,classified];
ArrayPlot[Sign[Abs[convolvedData]],DataRange->
{{Re[center]-1.2,Re[center]+1.2},{Im[center]-1.2,Im[center]+1.2}},
Epilog->{Red,Circle[{-0.12,0.13},0.07]}]

Pretty cool - but what’s with those do-hitchies, like the one circled in the image?

Big steps

The idea behind the big steps method is to compute some power series that approximate higher iterates of f, that is (Subscript[F, k])(z)≈(f^(2^k))(z) for z close to zero. More specifically, each (Subscript[F, k]) will be used on an annulus 1/(2^(k+K+1))< |z|<1/(2^(k+K)), where K is fixed. Typically, I find how many terms to use and on which annuli to use specific (Subscript[F, k])s experimentally.

Note that sometimes f already represents an iterate of the function you started with. In our harder set of examples, for example, we work the the function f(z)=(e^(2π i/n))z+(z^2), but (Subscript[F, 0]) is an approximation (f^n). Here’s some code to set this all up:

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

Here’s a little code to examine how well these polynomials work:

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

Looking at the fourth entry, for example, it appears that (Subscript[F, 4])(z) is a good approximation to (f^(7×(2^4)))(z) when z=1/(2^3) (and hopefully when |z|<1/(2^3)). Based on this, the following code defines a compiled function FF that uses (Subscript[F, k]) on the annulus 1/(2^k)<\[VerticalSeparator]z\[VerticalSeparator]≤1/(2^(k-1)) for k=4,5,…,10 and defaults to (f^n) for larger values of z.

FF = With[{
n = n,
f4=Function[z,Evaluate[F[4][z]]],
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]]]
},
Compile[{{z,_Complex}},
Which[
Abs[z]>1/2^3, 
Nest[Function[zz,N[Exp[2Pi*I/n]]zz+zz^2],z,n],
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]≤1/2^3,f4[z],
True,0]]];
iterate = With[{FF=FF,n=n},
Compile[{{z0,_Complex}},
Module[{z,i},
z=z0;
i = 0;
While[1/2^9<Abs[z]≤2 && i++<1000n,
z = FF[z]];
z],
RuntimeOptions -> "Speed", CompilationTarget -> "C"]];

We now compute some iteration data.

data = Monitor[
Table[iterate[x+I*y],{y,Im[center]+1.2,Im[center],-0.0025},
{x,Re[center]-1.2,Re[center]+1.2,0.0025}],y];//AbsoluteTiming

data = Join[data,Reverse[Rest[Reverse/@data]]];

And plot it as before.

kernel = {
{1,1,1},
{1,-8,1},
{1,1,1}
};
classified = Map[classify,data,{2}];
convolvedData = ListConvolve[kernel,classified];

ArrayPlot[Sign[Abs[convolvedData]]]

Definitely better!