Finding all points of period n in an iterated map

Question

I'm trying to implement an algorithm of Jenkinson and Pollicott to calculate the Hausdorff dimension of a Julia set for the map $f_c : z\mapsto z^2 + c$. It's described on page 40 of their paper, Calculating Hausdorff Dimension of Julia Sets and Kleinian Limit Sets. I can't seem to get the part about finding the periodic points with period $n$ to work.

The idea is to partition the line from $0$ to $c$ through the complex numbers and use Newton's method on each $c_i$ using the solution from one as the initial guess to the next. For $c=0$ we know the solution because they are roots of unity. We also take another precaution: to avoid the sensitivity to initial conditions that's prevalent in Newton's method, we find the roots of $F_{c,n} - I$ instead of finding the roots of $f^n_c - \operatorname{id}$, where $$F_{c,n}(z) = \langle f_c(z_n),f_c(z_1),\dots,f_c(z_{n-1})\rangle$$

I threw this together:

partitionSize = 10;
f[z_, c_] := z^2 + c;
fDim[z_, c_] := f[#, c] & /@ RotateLeft[z, -1];
findPeriodN[N_, c_] := 
  Nest[{FindRoot[fDim[z, #[[1]]] - z, {z, #[[2]]}] &, #[[2]] + 
      c/partitionSize &}, {NestList[#^2 &, (-1)^(1/(2 N - 1)), N], 0},
    partitionSize];
findPeriodN[3, 0.25];

And it doesn't work at all. Does anybody know a simple way to do this, or some other implementation that does the same thing?

Answer 1

Here's a technique that follows the paper's ideas fairly closely. We first implement f and F for orbit length 7 as follows.

Clear[f, F];
f[c_][z_] := z^2 + c;
F[c_][zs_List] := f[c] /@ RotateRight[zs];

c = I/4;
orbitOrder = 7;
vars = Table[z[k], {k, 1, orbitOrder}];

We can then use NSolve to find all orbits of this size straight away.

orbits = NSolve[F[c][vars] == vars, vars];
trimmed = Union[
   Union[#, SameTest -> (Abs[#1 - #2] < 10^(-12) &)] & /@ 
    Map[Last, orbits, {2}],
   SameTest -> (Length[#1] == Length[#2] && Norm[#1 - #2] < 10^(-8) &)];
trimmed // Length
trimmed[[1 ;; 5]] // Column

We've found 20 distinct orbits with the first two being actual fixed points. Let's check that the third is in fact periodic

z0 = trimmed[[3, 1]];
Nest[f[c], z0, 7] - z0 // Chop

(* Out: 0 *)

Now, the reference wants orbits of length up to 20 and this simple technique won't quite work, so let's try their root finding technique for length 21. The idea, as sketched in the paper, is to find periodic orbits for a sequence of c values $0=c_0,c_1,\dots,c_n$, where each step size is small and periodic orbits for $c_0$ are known. Thus, the periodic orbit for $c_i$ can be used as the initial guess for $c_{i+1}$. We can implement this as follows.

c = I/4;
orbitOrder = 21;
vars = Table[z[k], {k, 1, orbitOrder}];
init = NestList[#^2 &, Exp[1000*2 Pi*I/(2^orbitOrder - 1)], 
   orbitOrder - 1];
step = c/10;
lastGuess = init;
Do[
  newGuess = vars /.        
    FindRoot[F[i*step][vars] == vars, Transpose[{vars, lastGuess}]];
  lastGuess = newGuess,
  {i, 1, 10}];
orbit = NestList[f[c], newGuess[[1]], orbitOrder];
orbit // Column

The difficulty now is that this finds just one orbit at a time. Note, however, that the 1000 in the Exp is a parameter that affects the initial periodic orbit; thus, different choices can lead to differnt orbits for $c$.

Finally, here's a visualization of how the points jump around the Julia set:

waPic = WolframAlpha["julia set z^2 + i/4", {{"Plot", 1}, "Content"}];
pic = Graphics[{PointSize[0.001], Cases[waPic, _PointBox, Infinity] /. 
  PointBox -> Point}];
pts = {Re[#], Im[#]} & /@ orbit;

Show[{pic, Graphics[{Blue, {Arrowheads[0.05 Sqrt[Norm[First[#] - Last[#]]]], 
  Arrow[#]} & /@ Partition[pts, 2, 1], Red, PointSize[Medium], Point[pts]}]
}]

Answer 2

I've not seen the article but what you write suggests this is a crude homotopy from constant=0 to constant=c. This can be done more directly by setting up and solving a system of ODEs. Here is the general idea. We let constant vary from 0 to 'c' by multiplying by a time parameter t, force certain differential equations by differentiating with respect to t and setting that to zero. As noted in the query, we have initial conditions involving powers of a root of unity.

f[z_, c_] := z^2 + c;
fDim[z_, c_] := f[#, c] & /@ RotateLeft[z, -1];

findPeriodN[n_, c_] := Module[
  {t, z, vars, initval, sys, odes, ics, soln},
  vars = Map[#[t] &, Array[z, n]];
  initval = (-1)^(2/(2^n - 1));
  sys = fDim[vars, t*c] - vars;
  odes = Thread[D[sys, t] == 0];
  ics = Thread[(vars /. t -> 0) == NestList[#^2 &, initval, n - 1]];
  soln = NDSolve[Join[odes, ics], vars, {t, 0, 1}];
  vars /. soln /. t -> 1
  ]

Here is a simple run, for c=2 and we want a value with period 5.

res = findPeriodN[5, 2]

(* Out[234]= {{0.525826663857 + 1.29571894534 I, 
  0.597606069678 + 1.36264715528 I, 
  0.500325745942 + 1.6286524308 I, -0.402182908378 + 
   1.62971349102 I, -0.494214997214 - 1.31088584935 I}} *)

We check that the first one has period 5.

NestList[#^2 + 2 &, res[[1, 1]], 5]

(* Out[232]= {0.525826663857 + 1.29571894534 I, 
 0.597606095104 + 1.36264714065 I, 
 0.500325814984 + 1.62865247346 I, -0.402182958163 + 
  1.62971375222 I, -0.49421578233 - 1.31088619565 I, 
 0.525826621556 + 1.29572129346 I} *)

From here one could use FindRoot to polish that value, that is, improve it to a higher precision value, if so desired.