Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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},
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?

share|improve this question
This won't debug your code, but: (1) Use n rather than N. (2) I believe your starting point should be (-1)^(1/(2^n-1)) rather than (-1)^(1/(2*n-1)). – Daniel Lichtblau Feb 12 '13 at 15:38
Actually (-1)^(2/(2^n-1)). I'll post something with this in actual use. – Daniel Lichtblau Feb 12 '13 at 16:00
up vote 9 down vote accepted

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

enter image description here

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

enter image description here

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

enter image description here

share|improve this answer

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.

share|improve this answer
Very nice! The technique in the paper requires that we find all orbits up to period on the order of 20 or even more. The NDSolve command chokes on that but should run with the option Method->{"EquationSimplification"->"Solve"}. Also, I suppose that your initval could be fiddled with to allow for the need to find more orbits. – Mark McClure Feb 12 '13 at 16:21
@Mark McClure Also NDSolve might fare better with Method->Projection using the algebraic equations as the projection invariants. – Daniel Lichtblau Feb 12 '13 at 16:24

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.