Brown Method

Robert P. Munafo, 2003 Sep 22.

Here is a description of the algebraic methods used by A. Brown in 1981, John Stephenson in 1991-1992, Wolf Jung in 1996-1997, and others to derive exact formulas giving the locations of superstable points, bond points, and the boundaries of mu-atoms for periods 1 through 10.

As related by Wolf Jung, E. Netto published work similar to this in 1900. However, he was working with a different recurrence relation. In 1981 A. Brown was the first to derive the formulas for the Z²+c function, although he was working only in the reals. The work was substantially redone and generalized to the full complex case by Stephenson in 1991 and 1992. The algorithm has been refined and related to Galois theory (the same approach originally used by Netto) by Jung in 1996 and 1997.

The method creates polynomial equations relating the parameter value c to points in or on the unit disk D, mapping the disk to one or more mu-atoms of period n.

Technique

We start with the definition of the iterates, also called the recurrence relation:

Z_i+1 = Z_i² + c

where c is some parameter value whose iterates have period n. Z_i are the iterates of the the limit cycle of c, that is, Z₀ is the pre-periodic point and Z₀ = Z_n. For each c there is at most one Z₀ that satisfies this condition.

The goal is to find a polynomial g that relates D to c:

g(D, c) = 0

where D is a point in the unit disk. If such a polynomial is found, it can be used (either by explicit solving, as in the quadratic formula, or by numerical methods e.g. Newton's method) to locate the values of c for any given D or vice-versa.

This is possible because of a special property of the Mandelbrot set having to do with the derivatives of the iteration function: a point c is in the interior of a period n mu-atom, if and only if 2ⁿ times the n^th derivative of the n^th iteration Z_n has a modulus less than 1. This is a very fundamental property of the quadratic map, and is related to the connectedness of the Julia sets for stable points, to the way points in the interior of such sets converge asymptotically to a limit cycle, bifurcation, and so on.

We express this relationship with the equations:

P = Π_(1<=i<=n) [ Z_i ]

D = 2ⁿ P

P is the product of the n iterates and D is a point on the unit disk. To remove the factor of 2ⁿ, we substitute D = 2ⁿ P, and instead of g we seek a polynomial f in P and c:

f(P, c) = 0

As will be shown, this can be conducted with various algebra and abstract algebra techniques for each specific small period; so far it has been done for periods up to 10 (by Stephenson).

For periods above 2, Stephenson defines one or more of the following:

P_a = Π_(1<=i<=a) [ Z_i ] (for 0<a<n)

Additional special definitions are made later for even higher periods, as discussed below.

Throughout this entry, calculations are shown being performed with Maxima¹. Note that the input and output variables are all shown as %i and %o, omitting the numbers. I do not use them in the expressions although I do occasionally use % to refer to the recentmost output.

Period 1

For period 1, the derivation of f comes directly out of the recurrence relation. Here is the first iterate of a given point c:

(%i) z1=z0^2+c; 2 (%o) z1 = z0 + c

Since it's period 1, and we have stipulated that Z₀ is the pre-periodic point, it is equal to Z₁:

(%i) subst(z1,z0,%); 2 (%o) z1 = z1 + c

We want a polynomial in P and c, not in Z₁ and c. By the definition of P above, for period 1 P is equal to Z₁ so we can substitute directly:

(%i) subst(P,z1,%); 2 (%o) P = P + c

To get this in the canonical f(P,c) form, we shift things around like so:

(%i) reverse(%); 2 (%o) P + c = P (%i) % - P; 2 (%o) P - P + c = 0

So our polynomial f is f(P,c) = P² - P + c. We then substitute P=D/2 and solve for D and for c:

(%i) d6:subst(D/2,P,%); 2 D D (%o) -- - - + c = 0 4 2 (%i) solve(%, D); (%o) [D = 1 - SQRT(1 - 4 c), D = SQRT(1 - 4 c) + 1] (%i) solve(d6, c); 2 D - 2 D (%o) [c = - --------] 4 (%i) factor(%); (D - 2) D (%o) [c = - ---------] 4

Period 2

At period 2, we will start with the iterates and definition of P as before, adding the partial-product variable P₁ because there are now two Z's:

Z₂ = Z₁² + c

Z₁ = Z₂² + c

P = Z₁ Z₂

P₁ = Z₁

Here it is in Maxima:

(%i) E1:z2=z1^2+c; 2 (%o) z2 = z1 + c (%i) E2:z1=z2^2+c; 2 (%o) z1 = z2 + c

We substitute P₁ for Z₁ and P/P₁ for Z₂, giving two equations in 3 variables (P, P₁, and c) then eliminate the variable P₁ to get one equation in the two variables P and c:

(%i) E1:subst(P1,z1,subst(P/P1,z2,E1)); P 2 (%o) -- = P1 + c P1 (%i) E2:subst(P1,z1,subst(P/P1,z2,E2)); 2 P (%o) P1 = --- + c 2 P1 (%i) E1:ratsimp(P1 * E1); 3 (%o) P = P1 + c P1 (%i) E2:ratsimp(P1^2 * E2); 3 2 2 (%o) P1 = c P1 + P (%i) eliminate([E1,E2],[P1]); 2 4 3 2 2 2 4 3 2 (%o) [P (P - 3 P + (3 - 2 c ) P + (c - 1) P + c + 2 c + c )]

Note: eliminate uses resultants to determine how to algebraically manipulate a set of n+m input equations to produce n result equations that have m variables removed. If you want me to show how it is done explicitly, let me know by email.

(%i) ratsimp(first(%)); 6 5 2 4 2 3 4 3 2 2 (%o) P - 3 P + (3 - 2 c ) P + (c - 1) P + (c + 2 c + c ) P (%i) factor(%); 2 2 2 2 (%o) P (P - c - 1) (P + 2 c P - P + c )

Of the three terms shown here, the second is the one that actually corresponds to the period-2 mu-atom. (If this were not known, one could check the other term P² + 2 c P - P + c² numerically and one would find that it does not map the unit disk onto the period-2 mu-atom).

So we extract the term of interest, then convert from P back to D, and solve for c and for D:

(%i) first(first(rest(%))); (%o) P - c - 1 (%i) d11:subst(D/4,P,d10); D (%o) - - c - 1 4 (%i) solve(%,D); (%o) [D = 4 c + 4] (%i) solve(d11,c); D - 4 (%o) [c = -----] 4 (%i) expand(%); D (%o) [c = - - 1] 4

This very simple formula c=D/4-1 shows that the period-2 mu-atom R2.1/2a is a perfect circle of radius 1/4 and centered at -1 on the real axis.

Period 3

Here is the derivation for period 3. This time I am not bothering to go through the intermediate substitution of the P_i variables for the Z_i's — the algorithm used by Maxima's eliminate function works just as well on the Z_i equations as it would for P_i equations:

(%i) E1:z1=z3^2+c; 2 (%o) z1 = z3 + c (C2) E2:z2=z1^2+c; 2 (D2) z2 = z1 + c (C3) E3:z3=z2^2+c; 2 (D3) z3 = z2 + c (C4) E4:P=z1*z2*z3; (D4) P = z1 z2 z3 (C5) eliminate([E1,E2,E3,E4],[z1,z2,z3]); 6 8 7 3 6 3 5 (D5) [- P (P - 7 P + (4 c + 21) P + (- 7 c - 35) P 6 5 4 3 4 + (6 c + 6 c - 18 c - 11 c + 35) P 6 5 4 3 3 + (- c + 15 c + 51 c + 34 c - 21) P 9 8 6 5 4 3 2 + (4 c + 12 c - 28 c - 39 c - 45 c - 26 c + 7) P 9 8 6 5 4 3 12 11 10 + (- c - 3 c + 7 c + 9 c + 9 c + 5 c - 1) P + c + 6 c + 15 c 9 8 7 6 5 4 3 + 23 c + 27 c + 24 c + 16 c + 9 c + 3 c + c )] (C6) factor(first(%)); 6 2 3 2 3 2 3 (D6) - P (P - c P - 2 P + c + 2 c + c + 1) (P + 3 c P - P + c ) (C7) per3:first(first(rest(first(%)))); 2 3 2 (D7) P - c P - 2 P + c + 2 c + c + 1 (C8) per3:ratsimp(subst(a/8,P,per3)); 3 2 2 64 c + 128 c + (64 - 8 a) c + a - 16 a + 64 (D8) ---------------------------------------------- 64

Here, I have had to substitute a for P because Maxima groups the terms based on the alphabetical order of the variables and I want it to treat c as the most important variable. The form in which the solution is normally expressed looks like this:

c³ + 2 c² + (1-P) c + (1-P)² = 0 where P=a/8

In any case, the next step is to solve for either variable in terms of the other. Since the formula is of order 2 in P, the solution for P in terms of c comes from the quadratic formula; and since it is order 3 in c, the solution for c in terms of P comes from the cubic formula.

For explicit values, one can substitute and solve. For example, to find the bond points for R2.1/3.1/4a, R2.2/3.1/4a and R2F(1/2B1)S.1/4a we set a=0+i, because 0+i is the point which is 1/4 of the way around the unit disk, and we solve the formula as follows:

(C9) f:subst(%i,a,per3); 3 2 64 c + 128 c + (64 - 8 %I) c - 16 %I + 63 (D9) ------------------------------------------- 64 (C10) float(expand(float(rectform(solve(f=0,c))))); (D10) [c = 0.00977578406512 %I - 1.750480611554179, c = - 0.75087860018072 %I - 0.21919297274707, c = 0.7411028161156 %I - 0.03032641569875]

To get the general formulas to map points on the unit disk to points on R2.1/3a, R2.2/3a and R2F(1/2B1)Sa we substitute and solve as follows:

(C11) subst(D,a,expand(solve(per3,a))); (D11) [D = - 4 SQRT(- 4 c - 7) c + 4 c + 8, D = 4 SQRT(- 4 c - 7) c + 4 c + 8] (C12) subst(D,a,solve(per3,c)); SQRT(3) %I 1 (D10) [c = (- ---------- - -) 2 2 2 2 (D - 8) SQRT(27 D - 176 D + 1472) 27 D - 288 D + 1600 1/3 (---------------------------------- - --------------------) 384 SQRT(3) 3456 SQRT(3) %I 1 (---------- - -) (3 D + 8) 2 2 2 + ----------------------------------------------------------------- - -, 2 2 3 (D - 8) SQRT(27 D - 176 D + 1472) 27 D - 288 D + 1600 1/3 72 (---------------------------------- - --------------------) 384 SQRT(3) 3456 SQRT(3) %I 1 c = (---------- - -) 2 2 2 2 (D - 8) SQRT(27 D - 176 D + 1472) 27 D - 288 D + 1600 1/3 (---------------------------------- - --------------------) 384 SQRT(3) 3456 SQRT(3) %I 1 (- ---------- - -) (3 D + 8) 2 2 2 + ----------------------------------------------------------------- - -, 2 2 3 (D - 8) SQRT(27 D - 176 D + 1472) 27 D - 288 D + 1600 1/3 72 (---------------------------------- - --------------------) 384 SQRT(3) 3456 2 2 (D - 8) SQRT(27 D - 176 D + 1472) 27 D - 288 D + 1600 1/3 c = (---------------------------------- - --------------------) 384 SQRT(3) 3456 3 D + 8 2 + ----------------------------------------------------------------- - -] 2 2 3 (D - 8) SQRT(27 D - 176 D + 1472) 27 D - 288 D + 1600 1/3 72 (---------------------------------- - --------------------) 384 SQRT(3) 3456

There are three equations; the first two are for R2.1/3a and R2.2/3a and the third is for R2F(1/2B1)Sa

Period 4

For period 4, the technique that worked for periods 2 and 3 does not work. According to Stephenson (in Stephen1992A) this is because of "symmetries" in the equations (but it is not related to the fact that period 4 is a composite number and a multiple of period 2, because Stephenson reports similar problems for period 5). Essentially, the p+1 equations in p+2 variables do not actually contain p+1 distinct degrees of freedom of information, which they would need to contain in order to eliminate p variables and still have an equation left over. So, to eliminate the problem, Stephenson uses more equations.

Instead of the p Z_i, or the p P_i defined above, Stephenson starts with a family of functions S_m, for each m from 1 to the cycle period p, defined as follows:

S_m = SUM_1<=i<=p [ Z_i^m ]

or, specifically:

S₁ = Z₁ + Z₂ + Z₃ + ... + Z_p

S₂ = Z₁² + Z₂² + Z₃² + ... + Z_p²

...

S_p = Z₁^p + Z₂^p + Z₃^p + ... + Z_p^p

also, S₀ and all negative subscripts are defined to be 0. He also defines P_m by:

P_m = SUM{all combinations} [ Z_iZ_jZ_k ]

for example, when the period is 4, there are 4 P_m defined as follows:

P₁ = Z₁ + Z₂ + Z₃ + Z₄

P₂ = Z₁ Z₂ + Z₁ Z₃ + Z₁ Z₄ + Z₂ Z₃ + Z₂ Z₄ + Z₃ Z₄

P₃ = Z₁ Z₂ Z₃ + Z₁ Z₂ Z₄ + Z₁ Z₃ Z₄ + Z₂ Z₃ Z₄

P₄ = Z₁ Z₂ Z₃ Z₄

The formula for P₂ is the sum of all combinations of 2 Z_i's multiplied together, and the formula for P₃ is the sum of all combinations of 3 Z_i's multiplied together.

Jung's method yields the following relation between c and D:

6 5 4 3 2 c + 3 c + (1/16 D + 3) c + (1/16 D + 3) c - (1/16 D + 2) (1/16 D - 1) c 3 - (1/16 D - 1) = 0

Period 5

Jung's method yields the following:

15 14 13 12 11 c + 8 c + 28 c + (1/32 D + 60) c + (7/32 D + 94) c 2 10 / 11 2 33 \ 9 + (3/1024 D + 5/8 D + 116) c + |---- D + -- D + 114| c \1024 32 / 2 8 / 3 2 37 \ 7 + (3/512 D + 5/4 D + 94) c + |1/16384 D - 5/256 D + -- D + 69| c \ 32 / 2 6 + (3/32 D - 11) (3/1024 D - 3/32 D - 4) c / 3 2 33 \ 5 + (1/32 D - 1) |3/32768 D + 5/256 D - -- D - 26| c \ 32 / / 2 27 \ 2 4 + |3/1024 D + -- D + 14| (1/32 D - 1) c \ 32 / 3 3 4 2 - (3/16 D + 5) (1/32 D - 1) c + (1/32 D + 2) (1/32 D - 1) c 5 6 - (1/32 D - 1) c + (1/32 D - 1) = 0

Incomplete implementation in Maxima follows

/* Here I am trying to figure out how Pi and Si are defined as used in Stephen1992A */ /* First interpretation of Pi and Si */ z2:z1^2+c; p1:z1; p2:z1*z2; s1:z1; s2:p1^2-2*p2; z2:z1^2+c; p1:z1; p2:z1*z2; s1:z1+z2; s2:z1^2+(z1-c); p:2; z2:z1^2+c; p1:z1 + z2; p2:expand(z1*z2); s1:z1+z2; s2:z1^2+(z1-c); is(s2=p1 - p*c); p:3; z2:z1^2+c; z3:expand(z2^2+c); z4:expand(z3^2+c); reduce_z1(p):=expand(remainder(expand(p),z4)+z1*quotient(expand(p),z4)); p1:expand(z1+z2+z3); p2:expand(z1*z2+z3); p3:expand(z1*z2*z3); s1:expand(z1+z2+z3); /* s2:expand(z1^2+z2^2+(z1-c)); */ s2:reduce_z1(z1^2+z2^2+z3^2); /* s3:expand(z1^3+z2^3+z3*(z1-c)); */ s3:reduce_z1(z1^3+z2^3+z3^3); is(s2=p1 - p*c); q12:expand(s3+c*p1); p:4; z2:z1^2+c; z3:expand(z2^2+c); z4:expand(z3^2+c); z5:expand(z4^2+c); reduce_z1(p):=expand(remainder(expand(p),z5)+z1*quotient(expand(p),z5)); reduce_z2(p):=reduce_z1(expand(remainder(expand(p),z5^2)+z1^2*quotient(expand(p),z5*z5))); p1:expand(z1+z2+z3+z4); p2:expand(z1*z2+z3*z4); p2:reduce_z1(z1*z2+z2*z3+z3*z4+z4*z1); p3:expand(z1*z2*z3+z4); p4:expand(z1*z2*z3*z4); s1:expand(z1+z2+z3+z4); s2:expand(z1^2+z2^2+z3^2+(z1-c)); s2:reduce_z1(z1^2+z2^2+z3^2+z4^2); s3:expand(z1^3+z2^3+z3^3+z4*(z1-c)); s3:reduce_z1(z1^3+z2^3+z3^3+z4^3); s4:expand(z1^4+z2^4+(z1-c-2*c*z4+c^2)+(z1-c)^2); s4:reduce_z2(z1^4+z2^4+z3^4+z4^4); is(s2=expand(p1 - p*c)); is(s4=expand(p1*(1-2*c)+p*(c^2-c))); /* Now going back to lower periods to try to reverse-engineer the definition of the Pi: */ p:2; z2:z1^2+c; z3:expand(z2^2+c); reduce_z1(p):=expand(remainder(expand(p),z3)+z1*quotient(expand(p),z3)); p1:z1 + z2; p2:expand(z1*z2); s1:z1+z2; s2:z1^2+(z1-c); is(s2=p1 - p*c); p2rev:reduce_z1((p1^2-s2)/2); is(s2=reduce_z1(p1^2-2*p2)); /* Figuring out Q12 */ p:3; z2:z1^2+c; z3:expand(z2^2+c); z4:expand(z3^2+c); reduce_z1(p):=expand(remainder(expand(p),z4)+z1*quotient(expand(p),z4)); p1:expand(z1+z2+z3); /* p2:expand(z1*z2+z3); WRONG */ p2:reduce_z1(z1*z2+z2*z3+z3*z1); p3:expand(z1*z2*z3); s1:expand(z1+z2+z3); /* s2:expand(z1^2+z2^2+(z1-c)); */ s2:reduce_z1(z1^2+z2^2+z3^2); /* s3:expand(z1^3+z2^3+z3*(z1-c)); */ s3:reduce_z1(z1^3+z2^3+z3^3); is(s2=p1 - p*c); q12:expand(s3+c*p1); is(s2=reduce_z1(p1^2-2*p2)); is(ratsimp(q12)=ratsimp(z1*z2+z2*z3+z3*z1)); is(q12=expand(z1*z2+z2*z3+z3*z1)); is(s3=reduce_z1(p1^3-3*p2*p1+3*p3)); /* Now let's figure out P2 and P3 */ p:4; z2:z1^2+c; z3:expand(z2^2+c); z4:expand(z3^2+c); z5:expand(z4^2+c); reduce_z1(p):=expand(remainder(expand(p),z5)+z1*quotient(expand(p),z5)); reduce_z2(p):=reduce_z1(expand(remainder(expand(p),z5^2)+z1^2*quotient(expand(p),z5*z5))); p1:expand(z1+z2+z3+z4); p2:reduce_z1(z1*z2+z1*z3+z1*z4+z2*z3+z2*z4+z3*z4); p3:reduce_z1(z1*z2*z3+z1*z2*z4+z1*z3*z4+z2*z3*z4); p4:reduce_z1(z1*z2*z3*z4); s1:expand(z1+z2+z3+z4); s2:reduce_z1(z1^2+z2^2+z3^2+z4^2); s3:reduce_z1(z1^3+z2^3+z3^3+z4^3); s4:reduce_z2(z1^4+z2^4+z3^4+z4^4); is(s2=expand(p1 - p*c)); is(s4=expand(p1*(1-2*c)+p*(c^2-c))); q12:expand(s3+c*p1); p2rev:reduce_z2((p1^2-s2)/2); p2rev-reduce_z2(z3*z4)-reduce_z2(z2*z4)-reduce_z2(z1*z4) -reduce_z2(z2*z3)-reduce_z2(z1*z3)-reduce_z2(z1*z2); is(s2=reduce_z2(p1^2-2*p2)); is(ratsimp(q12)=ratsimp(z1*z2+z2*z3+z3*z4+z4*z1)); is(q12=expand(z1*z2+z2*z3+z3*z4+z4*z1)); is(s3=reduce_z1(p1^3-3*p2*p1+3*p3)); n:4$ print(" n:", n)$ (mo[1]:1, mo[2]:-1, mo[3]:-1, mo[4]:0)$ f[0]:z$ for i:1 thru n do ( f[i]:expand(f[i-1]^2+c), print(" f[", i, "] =", f[i]) )$ g:1$ for i:1 thru n do ( if (remainder(n, i)=0) then ( g:ratsimp(g * (f[i]-z)^(mo[n/i])), print(" g: ", g) ) )$ kn:hipow(g,z)/n$ print(" kn:", kn)$ q:z-a$ for i:1 thru n-1 do ( q:ratsimp(q+f[i]), print(" q:", q) )$ p:g$ print(" p: ", p)$ /* Might not always work -- might have to emulate "factors" instead */ while(hipow(q,z)>n) do ( r0:remainder(p, q, z), p:q, if(hipow(r0,z)>0) then ( q:r0, print(" q:", q) ) )$ r0:factor(num(remainder(p,q,z)))$ lr0:length(r0)$ for i:1 thru lr0 do ( term:first(r0), r0:rest(r0), if (hipow(term,z)=0) and (hipow(term,a)=kn) then ( h:term, print(" h:", h) ) )$

/* My attempt at translating the Maple code in Jung1997Some..Rational */ mobius(n) := block([f,rv,c], if (not integerp(n)) then return(FALSE), if (n=1) then return(1), if (primep(n)) then return(-1), rv:1, c:1, for f:2 thru n step 1 do ( if (primep(f)) then ( if (remainder(n,f)=0) then ( rv:-rv, n:n/f ), if (remainder(n,f)=0) then ( rv:0, n:n/f ), if (f > n) then return(rv) ) ), return(rv) ); for i:1 thru 20 step 1 do ( print("mobius[", i, "]:", mobius(i)) ); x(o) := block( if(atom(o)) then return([o]) else return(maplist('x,o)) )$ x(factor(12)); e1:expand((x^2+3*x+2)^5); x(e1); fac1(o) := block([t], t:factor(o), if(atom(t)) then return([t,1]) else return(maplist('fac1,t)) )$ factorlist(n) := (fac1(factor(n)))$ fac2b(i) := (i)$ flterm(o) := block([t], t:factor(o), if(atom(t)) then return([t,1]) else return(maplist('fac2b,t)) )$ factorlist(n) := block( if (n<0) then return(FALSE), if (n<2) then return([n,1]), if (primep(n)) then return([n,1]), maplist('flterm,factor(n)) )$ factorlist(0); factorlist(1); factorlist(2); factorlist(4); factorlist(12); factorlist(360); loprime(n) := block( if(not integerp(n)) then return(n), if(n<0) then return(loprime(-n)), if(n<4) then return(n), if(primep(n)) then return(n), return(loprime(first(factor(n)))) )$ for i:1 thru 20 step 1 do ( print("loprime[", i, "]:", loprime(i)) ); n:3$ print(" n:", n)$ f[0]:z$ for i:1 thru n do ( f[i]:expand(f[i-1]^2+c), print(" f[", i, "] =", f[i]) )$ g:1$ for i:1 thru n do ( print("loop2 ",i), if (remainder(n, i)=0) then ( g:ratsimp(g * (f[i]-z)^(mobius(n/i))), print(" g: ", g) ) )$ kn:hipow(g,z)/n$ print(" kn:", kn)$ q:z-a$ for i:1 thru n-1 do ( q:ratsimp(q+f[i]), print(" q:", q) )$ p:g$ print(" p: ", p)$ /* Might not always work -- might have to emulate "factors" instead */ while(hipow(q,z)>n) do ( r0:remainder(p, q, z), p:q, if(hipow(r0,z)>0) then ( q:r0, print(" q:", q) ) )$ r0:factor(num(remainder(p,q,z)))$ lr0:length(r0)$ for i:1 thru lr0 do ( print("loop4 ",i), term:first(r0), r0:rest(r0), if (hipow(term,z)=0) and (hipow(term,a)=kn) then ( h:term, print(" h:", h) ) )$

Footnotes

1 : Maxima, free software. http://maxima.sourceforge.net/

Mu-ency main page — index — recent changes — DEMZ

This page was written in the "embarrassingly readable" markup language RHTF, and was last updated on 2011 Jan 19.

s.27