mdbtxt1
mdbtxt2
Proceed to Safety

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 Z2+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:

Zi+1 = Zi2 + c

where c is some parameter value whose iterates have period n. Zi are the iterates of the the limit cycle of c, that is, Z0 is the pre-periodic point and Z0 = Zn. For each c there is at most one Z0 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 2n times the nth derivative of the nth iteration Zn 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) [ Zi ]

D = 2n P

P is the product of the n iterates and D is a point on the unit disk. To remove the factor of 2n, we substitute D = 2n 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:

Pa = Π(1<=i<=a) [ Zi ]                  (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 Maxima1. 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 Z0 is the pre-periodic point, it is equal to Z1:

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

We want a polynomial in P and c, not in Z1 and c. By the definition of P above, for period 1 P is equal to Z1 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) = P2 - 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 P1 because there are now two Z's:

Z2 = Z12 + c

Z1 = Z22 + c

P = Z1 Z2

P1 = Z1

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 P1 for Z1 and P/P1 for Z2, giving two equations in 3 variables (P, P1, and c) then eliminate the variable P1 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 P2 + 2 c P - P + c2 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 Pi variables for the Zi's — the algorithm used by Maxima's eliminate function works just as well on the Zi equations as it would for Pi 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:

c3 + 2 c2 + (1-P) c + (1-P)2 = 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 Zi, or the p Pi defined above, Stephenson starts with a family of functions Sm, for each m from 1 to the cycle period p, defined as follows:

Sm = SUM1<=i<=p [ Zim ]

or, specifically:

S1 = Z1 + Z2 + Z3 + ... + Zp

S2 = Z12 + Z22 + Z32 + ... + Zp2

...

Sp = Z1p + Z2p + Z3p + ... + Zpp

also, S0 and all negative subscripts are defined to be 0. He also defines Pm by:

Pm = SUM{all combinations} [ ZiZjZk ]

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

P1 = Z1 + Z2 + Z3 + Z4

P2 = Z1 Z2 + Z1 Z3 + Z1 Z4 + Z2 Z3 + Z2 Z4 + Z3 Z4

P3 = Z1 Z2 Z3 + Z1 Z2 Z4 + Z1 Z3 Z4 + Z2 Z3 Z4

P4 = Z1 Z2 Z3 Z4

The formula for P2 is the sum of all combinations of 2 Zi's multiplied together, and the formula for P3 is the sum of all combinations of 3 Zi'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/




From the Mandelbrot Set Glossary and Encyclopedia, by Robert Munafo, (c) 1987-2024.

Mu-ency main pageindexrecent changesDEMZ


Robert Munafo's home pages on AWS    © 1996-2024 Robert P. Munafo.    about    contact
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License. Details here.

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