Parabolic points

A neutral fixed point is a fixed point (Subscript[z, 0]) where |f'((Subscript[z, 0]))|=1. Investigating dynamic iteration near a neutral point is generally difficult because things move very slowly there.

A neutral fixed point is called parabolic if f'((Subscript[z, 0]))=(e^(α π i)), where α is a rational number. There is additional structure near a neutral fixed point that illuminates the behavior and helps us generate nice images.

A challenge

We’re going to begin this with a challenge - generate the Julia set of f(z)=z+(z^5)

f(z)=z+(z^5).

A simple polynomial - how hard could it be?? Here are a couple of images using escape time and inverse iteration algorithms.

GraphicsRow[{
JuliaSetPlot[z+z^5,z],
JuliaSetPlot[z+z^5,z, ColorFunction->None]
}, ImageSize->600]

It looks like there’s some weird hole near zero but it’s actually easy to prove that zero is in the Julia set.

The difficulty

Suppose that zero is a parabolic point for f with 0<z<f(z) for small, positive z.  
Starting at (Subscript[z, 0])=a.  How many iterates does it take to reach b?

Supposing it takes n times and letting (Subscript[z, i])=(f^i)((Subscript[z, 0])), then we can write

n=(∑_(i=0)^(n-1))(1/(f((Subscript[z, i]))-(Subscript[z, i])))((Subscript[z, i+1])-(Subscript[z, i]))≈(∫_a^b)(1/(f(z)-z))dz.

If we let a=1/64, b=1/32, and f(z)=z+(z^5), then we get

(∫_(1/64)^(1/32))1/(z^5)\[DifferentialD]z

The exact value obtained by running the loop is 3932162.

zi = 1/64.0;
cnt = 0;
While[zi≤1/32.0,
cnt=cnt+1;
zi= zi+zi^5];
cnt

A solution

Here’s a more detailed look at our naive approach.

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

It turns out that if we iterate from inside a bounded region of the Fatou set, then we always end up in one of the four lobes that meet at the origin. We can shade the initial seed by which lobe we end up in.

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]]]];
classified = Map[simpleClassify[#[[1]]]&,data,{2}];
classifiedPic = ArrayPlot[classified, DataRange -> {{-1,1},{-1,1}}]

This works because there are four weakly attractive directions and four weakly repulsive directions.

f[z_] = z+(z^5);
Show[classifiedPic, Epilog -> {
{Red,
Table[Arrow[{Re[#],Im[#]}& /@ 
{0.75Exp[k*Pi*I/4],f[0.75Exp[k*Pi*I/4]]}],
{k,0,6,2}]},
{Green,
Table[Arrow[{Re[#],Im[#]}& /@ 
{0.75Exp[k*Pi*I/4],f[0.75Exp[k*Pi*I/4]]}],
{k,1,7,2}]}
}]

You can see this by examining f(r (e^(k π i/4))):

f(r (e^(k π i/4))) = r (e^(k π i/4))+(r^5) (e^(5k π i/4)) = r (e^(k π i/4))+(r^5) (e^(k π i/4)), if k is even
f(r (e^(k π i/4))) = r (e^(k π i/4))+(r^5) (e^(5k π i/4)) = r (e^(k π i/4))-(r^5) (e^(k π i/4)), if k is odd

Thus, we have a little step away from the origin when r is small and k is even but we have a little step toward the origin when r is small and k is odd.

Now, the Julia set is the boundary between those regions so we use a technique called 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]]

The Leau-Fatou flower theorem

The theorem

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]

A more general classifier

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

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'/@ %]

Let’s figure out the nature of the parabolic point at zero.

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

That indicates three fold symmetry. We should iterate while we’re bounded by 2 away from the attractive fixed point 0.8+0.8i; we’ll return the landing value.

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

We’ll then classify on whether we diverged, landed near the attractive fixed point, or by which lobe we landed in.

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}];
ArrayPlot[classified]

And finally boundary scan.

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

A harder set of examples

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

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]

Not too bad.

Big steps

An even harder example

Now let’s try f(z)=(e^(2π i/7))z+(z^2).

n=21;
f[z_] = N[Exp[2Pi*I/n],20]z+z^2;
F[z_] = z;
Do[
F[z_] = Normal[Series[F[f[z]],{z,0,n+1}]],
{n}];
F[z_] = Chop[N[F[z]]];
Clear[v,vv];
a = Coefficient[F[z],z,n+1];
vv = v/.Solve[n*a*v^n==1,v];
v = First[vv]

series

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]

Arg[%]/Pi

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.005},
{x,Re[center]-1.2,Re[center]+1.2,0.005}],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?

The big steps algorithm

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,40}]]];
Do[F[0][z_] = Chop[N[Normal[Series[F[0][f[z]],{z,0,40}]]],10^-5],
{n-1}];
Do[F[k][z_] = Chop[N[Normal[Series[F[k-1][F[k-1][z]],{z,0,40}]]],10^-5],
{k,1,10}]

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],
CompilationOptions->{"InlineCompiledFunctions"->False},
RuntimeOptions -> "Speed",CompilationTarget -> "C"]];

iterate[0.01+0.02I]

0.009999999986557628` +0.019999999973115256` I

0.00999999999999936` +0.01999999999999872` I

iterate[[6]]

We now compute some iteration data.

(0.917622 + 0.147377)/0.005

%/60

data = Monitor[
Table[iterate[x+I*y],{y,Im[center]+1.2,Im[center],-0.005},
{x,Re[center]-1.2,Re[center]+1.2,0.005}],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]]]

NotebookSave[]

Definitely better!