# Mandelbrot boundary

Is there a sequence of parameterized expressions for the border of all the major bulbs of the mandelbrot set? By major meaning all bulbs with diameter greater than 0.01 for example. I am interested in generating the major features of the set to apply transformations to it. I want to create a vector of points of the border (for major features only) for a Matlab type language. I know the parameterized expression for the main cardiod and the period 2 bulb centered at -1.

In the discussion, the op asked wheather the bulbs are exactly circular, and if I could post the code to generate the Taylor series for the bulbs. The quick answer is that only the period2 bulb centered at -1 is exactly circular. My first answer gave the center and size of some of the bulbs and cardioids. The size is a rather crude approximation of a bulb or cardiod, since all but one of the bulbs isn't exactly circular, and since the cardioids also need a direction. I enhanced the original pari-gp code to also calculate the mapping of the bulbs/cardioid to the unit circle. The pari-gp code is below. islopetaylor(c,n) returns the taylor series for the period n cardiod or bulb near c.

All of these hyperbolic centers can have their bulbs and/or cardioids mapped to the unit circle, where the boundary of the unit circle corresponds to where the $|\text{slope}|=1$. The inside of the circle, is mapped to points inside the hyperbolic region, where the region has an attracting fixed point of period n.

Finding such a region involves getting an approximation for the nearest fixed point, in the neighborhood of the hyperbolic center. The fixed point is attracting, so we can iterate, starting with z=0 with period n, $z \mapsto f^{n}(z)$, where $f(z)=z^2+c$, and c is the point whose fixed point is desired, so that in the limit, the fixed point is $l=\lim_{k \to \infty} f^{kn}(0)$, $f^{n}(l)=l$.

At the fixed point l, calculate the derivitive. $\frac{dy}{dx}f^{n}(l+x)$, where l is the fixed point, described above. If we call this derivitive fixed point function g(c), then the desired function mapping the bulb/cariod to the unit circle is $g^{-1}(z)$. For a given bulb, $g^{-1}(0)$ is hyperbolic center of that bulb or cardiod. I generated the Taylor series by using a circular Cauchy integral numerical approximation at a radius of 0.5, using 50 sample points; which gives results accurate to approximately 20 decimal digits or so, with a 25 term Taylor series.

Here is an example of how to use the pari-gp code, below, to print out the taylor series for the period4 bulbs and cardioids. EDIT: I also included the results for all of the period3 and period5 bulbs and cardioids. Similar code would work for other periods, and sizes. The prtpoly routine is just for pretty printing, of the first few terms; The polynomial output of of islopetaylor can be used directly. For bulbs, the diameter is approximately twice the absolute value of the x coefficient; the Taylor series is obviously much more exact. My first answer estimate of 1.6x the distance to the nearest 2n bulb is also a good approximation. p=periodnzero(n,l); also works for n=6..10. There are a total of sixteen period6 bulbs and cardioids, most of which are pretty small. p=randomperiodn(n,l) works for the larger bulbs, for p>=5.

p=periodnzero(3,0.001); /* rtn period3 bulbs size>0.001 */
/* taylor series for p3 bulbs, first 7 terms of gt poly  */
for (n=1,length(p),gt=islopetaylor(p[n]); prtpoly(gt,7));
p=periodnzero(4,0.001); /* rtn period4 bulbs size>0.001 */
for (n=1,length(p),gt=islopetaylor(p[n]); prtpoly(gt,7));
p=periodnzero(5,0.001); /* rtn period5 bulbs size>0.001 */

---- below is the output ----

(08:49) gp > p=periodnzero(3,0.001); /* rtn period3 bulbs size>0.001 */
periodnzero(3);
0.0288243794158896630 -1.75487766624669276
0.185709338183664371 -0.122561166876653620 + 0.744861766619744237*I
(09:06) gp > /* taylor series for p3 bulbs, first 7 terms of gt poly  */
(09:06) gp > for (n=1,length(p),gt=islopetaylor(p[n]); prtpoly(gt,7));
{z=   -1.75487766624669276
+x^ 1* 0.00951775795656622345
+x^ 2*-0.00439215019285891483
+x^ 3*-0.000255593673326239667
+x^ 4* 0.00000507529911656475223
+x^ 5* 0.00000242235932732248524
+x^ 6* 0.000000167831763783128979
}
{z=   (-0.122561166876653620 + 0.744861766619744237*I)
+x^ 1*(-0.00475887897828311173 - 0.0943369626738123763*I)
+x^ 2*( 0.00219607509642945741 - 0.00111573891029371424*I)
+x^ 3*( 0.000127796836663119834 + 0.0000966580666527584112*I)
+x^ 4*(-0.00000253764955828237611 + 0.0000131076055547205131*I)
+x^ 5*(-0.00000121117966366124262 + 0.000000329634003020475535*I)
+x^ 6*(-0.0000000839158818915644933 - 0.0000000951268060412508001*I)
}
(09:06) gp > p=periodnzero(4,0.001); /* rtn period4 bulbs size>0.001 */
periodnzero(4);
0.00157165406194145265 -1.94079980652948475
0.113351748952365738 -1.31070264133683288
0.0129942224028632946 -0.156520166833755062 + 1.03224710892283180*I
0.0868986952447806464 0.282271390766913880 + 0.530060617578525299*I
(09:06) gp > for (n=1,length(p),gt=islopetaylor(p[n]); prtpoly(gt,7));
{z=   -1.94079980652948475
+x^ 1* 0.000495813881701186015
+x^ 2*-0.000244518035475789364
+x^ 3*-0.00000233355263604428247
+x^ 4* 0.0000000534459067501123460
+x^ 5* 0.00000000181828170546195257
+x^ 6*-3.23710018614877151 E-12
}
{z=   -1.31070264133683288
+x^ 1* 0.0589801138156112997
+x^ 2* 0.00165047843656841052
+x^ 3* 0.0000693045614883186867
+x^ 4* 0.00000269896930057060465
+x^ 5* 0.0000000523314280419260434
+x^ 6*-0.00000000565936605551606439
}
{z=   (-0.156520166833755062 + 1.03224710892283180*I)
+x^ 1*( 0.00345889762902884997 - 0.00245474033761904195*I)
+x^ 2*(-0.00152469430548367362 + 0.00130179544716022133*I)
+x^ 3*(-0.000143828841409685183 - 0.0000378568944213399944*I)
+x^ 4*( 0.00000385828834515733905 - 0.00000994102399364532939*I)
+x^ 5*( 0.00000140472771445005478 + 0.000000673772055503541835*I)
+x^ 6*(-0.0000000523614237035883056 + 0.000000196431825642151357*I)
}
{z=   ( 0.282271390766913880 + 0.530060617578525299*I)
+x^ 1*(-0.0331968614776850928 - 0.0287926429442936276*I)
+x^ 2*( 0.000821714104937363043 - 0.00133634537333004157*I)
+x^ 3*( 0.000110343336983547981 + 0.0000572780596940444853*I)
+x^ 4*(-0.00000523449594881769755 + 0.0000117988699179895988*I)
+x^ 5*(-0.00000143180256932374878 - 0.000000528764025646146377*I)
+x^ 6*( 0.0000000551927252814394115 - 0.000000186571532024463594*I)
}
(09:06) gp > p=periodnzero(5,0.001); /* rtn period5 bulbs size>0.001 */
periodnzero(5);
0.00124142647819312318 -1.86078252220485487
0.00643004051683055908 -1.62541372512330374
0.0782244856281390343 -0.504340175446244000 + 0.562765761452981964*I
0.00750234541167943880 0.359259224758007439 + 0.642513737138542349*I
0.0471208943401874649 0.379513588015923745 + 0.334932305597497587*I
0.00540957192966561334 -0.0442123577040706231 + 0.986580976280892768*I
0.00109467916151353002 -0.198042099364253840 + 1.10026953729269853*I
0.00454737524002394192 -1.25636793006818076 + 0.380320963472722507*I
(09:06) gp > for (n=1,length(p),gt=islopetaylor(p[n]); prtpoly(gt,7));
{z=   -1.86078252220485487
+x^ 1* 0.000390175993253921407
+x^ 2*-0.000193606895721835075
+x^ 3*-0.00000111641882446931317
+x^ 4* 0.0000000945861637002427640
+x^ 5* 0.00000000182079818330302844
+x^ 6*-6.26858107084452022 E-11
}
{z=   -1.62541372512330374
+x^ 1* 0.00202943297897757605
+x^ 2*-0.00100562087600704118
+x^ 3*-0.0000100969021567375794
+x^ 4* 0.00000299007153447035234
+x^ 5* 0.0000000448549883194611938
+x^ 6*-0.0000000140393194397544057
}
{z=   (-0.504340175446244000 + 0.562765761452981964*I)
+x^ 1*( 0.0227804277538978753 - 0.0313962947316761485*I)
+x^ 2*(-0.000187493753264388258 + 0.000240573004638874982*I)
+x^ 3*(-0.0000173170154681859400 + 0.0000483439272385087499*I)
+x^ 4*( 0.00000186140591725642208 - 0.00000140612260499651276*I)
+x^ 5*(-0.0000000359687742865028775 - 0.000000238921821822236246*I)
+x^ 6*(-0.0000000140479390745627075 + 0.0000000154228924213440731*I)
}
{z=   ( 0.359259224758007439 + 0.642513737138542349*I)
+x^ 1*( 0.00125062477418773157 - 0.00205512908763713047*I)
+x^ 2*(-0.000473500958540455897 + 0.00108015520081310997*I)
+x^ 3*(-0.000108144537904337611 - 0.0000221136751964296838*I)
+x^ 4*( 0.00000365352385186322974 - 0.0000109205152171821385*I)
+x^ 5*( 0.00000143644262907468698 + 0.000000753636636313957017*I)
+x^ 6*(-0.000000133837667274211728 + 0.000000211504575842695306*I)
}
{z=   ( 0.379513588015923745 + 0.334932305597497587*I)
+x^ 1*(-0.0230612323455722479 - 0.00538785818574664379*I)
+x^ 2*( 0.000223715667692254576 - 0.00100096226615781152*I)
+x^ 3*( 0.0000924390752585046801 + 0.0000283948327212593745*I)
+x^ 4*(-0.00000444777220263711731 + 0.0000109991366640167731*I)
+x^ 5*(-0.00000147488363872001793 - 0.000000750149091456313290*I)
+x^ 6*( 0.000000131362683170454543 - 0.000000211363105809565828*I)
}
{z=   (-0.0442123577040706231 + 0.986580976280892768*I)
+x^ 1*(-0.00168631410697906826 + 0.000419905880430706866*I)
+x^ 2*( 0.000824122737737012684 - 0.000147452933665481287*I)
+x^ 3*( 0.0000148294436540590271 - 0.0000423012304388835166*I)
+x^ 4*(-0.00000179611219530982926 + 0.000000252886482712731910*I)
+x^ 5*( 0.000000132428561354388846 + 0.000000194421657850058569*I)
+x^ 6*( 0.0000000186904817675855266 - 0.0000000121128281556966309*I)
}
{z=   (-0.198042099364253840 + 1.10026953729269853*I)
+x^ 1*( 0.000338411530274354574 + 0.0000551258574846720436*I)
+x^ 2*(-0.000169328304150344563 - 0.0000234867261464734632*I)
+x^ 3*( 0.0000000298451123227361170 - 0.00000280280755567260578*I)
+x^ 4*( 0.0000000425917226359921085 + 0.0000000627926935964822686*I)
+x^ 5*(-0.00000000297733729350064710 + 0.00000000108485294983731985*I)
+x^ 6*( 6.65684728185643910 E-12 - 9.88920256175429103 E-11*I)
}
{z=   (-1.25636793006818076 + 0.380320963472722507*I)
+x^ 1*(-0.000846269215461327949 - 0.00116600657512957892*I)
+x^ 2*( 0.000389349729034723316 + 0.000591786929496725821*I)
+x^ 3*( 0.0000237849154393727229 - 0.00000452582161657443604*I)
+x^ 4*(-0.000000856099479577055100 - 0.00000108478239958016467*I)
+x^ 5*(-0.0000000783802180386627187 + 0.0000000644312584363213421*I)
+x^ 6*( 0.00000000487678960666910970 + 0.00000000500884893452725251*I)
}
(09:06) gp > p=periodnzero(5,0.001); /* rtn period3 bulbs size>0.001 */
periodnzero(5);
0.00124142647819312318 -1.86078252220485487
0.00643004051683055908 -1.62541372512330374
0.0782244856281390343 -0.504340175446244000 + 0.562765761452981964*I
0.00750234541167943880 0.359259224758007439 + 0.642513737138542349*I
0.0471208943401874649 0.379513588015923745 + 0.334932305597497587*I
0.00540957192966561334 -0.0442123577040706231 + 0.986580976280892768*I
0.00109467916151353002 -0.198042099364253840 + 1.10026953729269853*I
0.00454737524002394192 -1.25636793006818076 + 0.380320963472722507*I


Here is the pari-gp program to generate this. This code also includes updates for the code in the original answer well. More updates, including the plotsetup routine, and maketheplot routine; see the picture below.

print ("periodnzero(n,l); /* calculate all periodn bulbs size>l, returns vector */");
print ("randomperiodn(n,l); /* n>4 approximate all periodn bulbs size>l rtn vect */");
print ("invzero(c,n); /* use Newton's method */");
print ("estim2nzero(c,n); /* find nearest 2n zero */");
print ("gt=islopetaylor(c,n); /* taylor series for c, n iterations */ ");
print ("prtpoly(gt,10); /* print taylor series for islopetaylor */ ");
default(format,"g0.18");
z=1.0;
precis=precision(z);
plim=10^(-0.5*precis); /* precision limit */
/* plim=10^(-0.65*precis); */ /* precision limit */
periodnzero(n,l) = {
local(z,i,pcur,otemp,Cp,oout);
otemp=vector(2^(n-1));
pcur=0;
zn=0;
print("periodnzero("n");");
Cp=croot(n);
for (i=1,length(Cp),
z=estimz2(Cp[i],n,l);
if (z<>0,pcur++;otemp[pcur]=Cp[i]);
);
oout=vector(pcur);
i=1;while (otemp[i]<>0,oout[i]=otemp[i];i++);
return(oout);
}
/* return all the imag>0 roots of zn, iterating zn<=zn^2+x */
croot(m) = {
local(zn,Cr,n,i,v);
zn=x;
for (n=1,m-1, zn=zn^2+x; );
Cr=polroots(zn);
i=0;
v=vector(length(Cr));
for (n=1,length(Cr),
z=Cr[n];
if (((real(z)<>0) || (imag(z)<>0)) && (imag(z)>=0),
if (imag(z)==0, z=real(z));
i++;
v[i]=z;
);
);
return(v);
}
estimz2(C0,n,l) = {
local(z,y,i);
z=0;
y=0;
if (l==0,l=0.01);
/* check to make sure that the z_n root is not repeating by a factor of n */
for (i=1,n-1,z=z^2+C0;if (abs(z)<plim,y=i));
if (y==0,
y=estim2nzero(C0,n);  /* estimate for 2n zero */
y=invzero(y,n*2); /* refinement for estimate */
/* 1.647(zn-z2n) is approx size of the bug or bulb, based on Feignenbaum limit */
if (1.6*abs(C0-y)>l,
print(1.6*abs(C0-y)" "C0 );
return(y);
);
return(0);
,
return(0);
);
}
sizez2(C0,n) = {
1.6*abs(C0-invzero(estim2nzero(C0,n),n*2));
}
/* use newton's method, centered at C0 */
/* iterating z=z^2+x+C0 with z only needing a0 and a1 terms */
/* s<>0 is used to calculate inv(s) instead of inv(0) */
/* m<>0 is used to calculate Misiurewicz point near C0,n,m */
invzero(C0,n,s,m)={
local(z,i,j,a1,sz,sa1);
z=1;
j=0;
while ( ((abs(z)>plim)||(j<1)) && (abs(z)<2) && (j<50),
j++;
z=0;  a1=0;
sz=0; sa1=0; /* Misiu */
i=0;
while ((i<n) && (abs(z)<2),
a1=2*a1*z+1;
z=z^2+C0;
i++;
if (i==m, sz=z; sa1=a1; ); /* Misiu */
);
z=z+sz-s;
a1=a1+sa1;
C0=C0-z/a1;
/*  if (j<3, C0=C0-z/(2*a1), C0=C0-z/a1); */
/*  C0=C0-(z-s)/a1; */
/*  if ((abs(z-s)>plim),print(z-s)); */
);
jm=0.001*a1^-2;
if (abs(z)>plim, return(0));
return(C0);
}
/* find nearest 2n zero, assuming that n is a zero */
estim2nzero(C0,n) = {
local(z1,z2,y,i,a0,a1,a2,a3);
if (C0==0,return(-1));
a0=0;a1=0;a2=0;a3=0;
/* iterating 2n times, centered at C0, */
/* starting with z=0; z<=z^2+x+C0, truncating to x^3 terms */
for (i=1,2*n,
a3=2*a3*a0+2*a1*a2;
a2=2*a2*a0+a1^2;
a1=2*a1*a0+1;
a0=a0^2+C0;
);
/* assume zn=0, so z2n=0, so a0=0, */
/* divide by x, solve quadratic root closest to zero */
a0=a1;
a1=a2;
a2=a3;
y=a1^2-4*a0*a2;
z1=(-a1+sqrt(y))/(2*a2);
z2=(-a1-sqrt(y))/(2*a2);
if (abs(z1)>abs(z2),z1=z2);
return(C0+z1);
}

shortperiodn(n,nc,cnt,l) = {
local(z,y,i,j);
z=0;for (i=1,n,z=truncs(z^2-0+x+nc,cnt+1));
v=polroots(z);
for (i=1,cnt,
z=v[i]+nc;
y=0;
for (j=1,n,y=y^2+z);
y=abs(y);
if (y<0.04,
z=invzero(z,n);
if (abs(imag(z))<1E-21, z=real(z));
if (imag(z)<0,z=conj(z));
if (((real(z)<>0) || (imag(z)<>0)) && (imag(z)>=0),
j=1;
while (abs(oj[j]-z)>1E-21 && (oj[j]<>0), j++);
if (oj[j]==0,
if (estimz2(z,n,l)<>0, oj[j]=z; if (j==15, oj[j]=v[i]+nc));
);
);
);
)
}
randomperiodn(n,l,j)= {
local(i,y,oout);
if (l==0,l=0.01);
if (j==0, j=50);
oj=vector(128);
print("randomperiodn("n");");
for (i=1,j,
shortperiodn(n,cfromz(goldr^i),16,l);    /* main cardiod */
shortperiodn(n,-1+0.25*(goldr^i),16,l);  /* 2n bulb */
shortperiodn(n,-2+fwidth*((gold*i)%1),16,l); /* tip */
);
i=1;while (oj[i]<>0,i++);
oout=vector(i-1);
i=1;while (oj[i]<>0,oout[i]=oj[i];i++);
return(oout);
}
/* truncate series */
truncs(z,n)={
local(y,i);
y=0;
if (n==0,n=8);
for (i=0,n-1,y=y+polcoeff(z,i)*x^i);
return(y);
}
gfunc(C0,n)={local(i,z);z=0;for (i=1,n,z=z^2+C0);z}
cfromz(z) = subst(x/2-x^2/4,x,z);
lfromc(C0) = {1/2-sqrt(1/4-C0);}
gold=2/(sqrt(5)+1);
goldr=conj(exp(2*Pi*I*gold));
/* feigenbaum point+2, and feigenbaum ratio */
fwidth= 0.598844810907949399476173212106138707773702;
feign = 4.669201609102990671853203820466201617258186;

/* slope at the fixed point of C0, using n iterations */
slopeh(C0,n)={
local(l,l0,l1,i,j,co,s);
s=1;
co=0;
while (abs(s)>plim, /* one past plim */
l0=co;l1=1; /* l1=slope */
for (i=1,n,l0=l0^2+C0;l1=l1*2*l0);
s=l0-co; /* as co approaches fixed point, (l0-co) approaches zero */
co=co-(l0-co)/(l1-1); /* co=updated approximation for fixed point */
);
l0=co;l1=1; /* one more iteration to improve l1 slope approximation */
for (i=1,n,l0=l0^2+C0;l1=l1*2*l0);
return(l1);
}
invslopeh(z,C0,n,est) = {
local (y,s,slop,lastyz,curyz,lest,ly,pgoal);
lastyz=100;
y=slopeh(C0+est,n);
curyz=abs(y-z);
lest=est+curyz*jm;
ly=slopeh(C0+lest,n);
pgoal=10^(-precis/1.3);
/* generate the fixed point for pentation by iteration slog */
s=1;
while ((curyz>pgoal) && ((curyz<lastyz) || (s<3)),
est=precision(est,precis);
y=precision(y,precis);
slop=(y-ly)/(est-lest);
lest=est;
ly=y;
est=est+(z-y)/slop;
lastyz=curyz;
y=slopeh(C0+est,n);
curyz=abs(y-z);
s++;
);
if (curyz>0.1, print (curyz " bad result, need better initial est"));
return(est);
}
/* default use r=0.5, samples=50 */
islopetaylor(C0,n,r,samples) = {
local(rinv,s,t,x1,y,z,tot,t_est,tcrc,halfsamples,wtaylor,terms);
if (n==0,
n=1;
z=C0;
while (abs(z)>plim,z=z^2+C0;n++);
);
if (r==0,r=0.5);
C0=invzero(C0,n);
if (samples==0, samples=50);
halfsamples=samples/2;
terms = floor(samples*0.51);
t_est    = vector (samples,i,0);
tcrc     = vector (samples,i,0);
if (r==0,r=1);
rinv = 1/r;
wtaylor=C0;
for(s=1, samples, x1=-1/(samples)+(s/halfsamples); tcrc[s]=exp(Pi*I*x1); );

for (t=1,samples, t_est[t] = invslopeh(r*tcrc[t],C0,n); );

for (s=0,terms-1,
tot=0;
for (t=1,samples,
tot=tot+t_est[t];
t_est[t]=t_est[t]*conj(tcrc[t]);
);
tot=tot/samples;
if (s>=1, tot=tot*(rinv)^s);
wtaylor=wtaylor+tot*x^s;
);
wtaylor=precision(wtaylor,precis);
if (imag(C0)==0, wtaylor=real(wtaylor));
return(wtaylor);
}
prtpoly(wtaylor,t) = {
local(s,z,iprt);
if (t==0,t=7);
z=polcoeff(wtaylor,0);
if (imag(z)<>0,iprt=1,iprt=0);
print1 ("{z=   ");
if (iprt,print1("("));
if (real(z)<0, print1(z), print1(" " z));
if (iprt,print(")"),print() );
for (s=1,t-1,
z=polcoeff(wtaylor,s);
if (s>9, print1("+x^" s), print1("+x^ " s));
if (iprt,print1("*("),print1("*") );
if (real(z)<0, print1(z), print1(" " z));
if (iprt,print(")"),print() );
);
print("}");
}

local(n,i,zeros,zeron);
zp=0;
if (numb==0,numb=85);
zeros=vector(numb);zeron=vector(numb);zerov=vector(numb);
if (low==0,low=1);
if (high==0, high=14);
for (n=low,high,
for (i=1,length(C),
zp++;
zeros[zp]=C[i];
zeron[zp]=n;
);
);
ip=0;
print(zp);
for (i=1,zp,
ip++;
zerov[ip]=truncs(islopetaylor(zeros[i],zeron[i]),7);
if (imag(zeros[i])<>0,
ip++;
zerov[ip]=conj(zerov[ip-1]);
);
);
}
evaln(z)={local(y,s,i);
y=vector(ip*2);
for (i=1,ip,
s=subst(zerov[i],x,z);
y[i*2-1]=real(s);
y[i*2]=imag(s);
);
return(y)
}
maketheplot(w) = local(t,z); { ploth(t=0,2*Pi,z=exp(t*I);evaln(z),1); }


Mandelbrot plot using bulbs/cardioids from period 1..20, bigger than 0.004. I truncated the Taylor series to seven terms, and plotted 215 different bulbs and cardioids; 112 are unique if not counting complex conjugates twice. For the randomperiodn, I used (n,0.004,350). It took 350 random tries to get both of the period(15) sub-bulbs on the period(5) bulb; some bulbs are still missing. Directly searching for bulbs near $\frac{k\pi i}{n}$ and other algorithms would probably work better; I'm pretty new at this.

• Thanks again. Are these algorithms based on theory from the book you recommended? – PMay Sep 22 '13 at 19:36
• The first answer used Devaney's article's description of how to calculate the hyperbolic center of the bulbs and cardiods. The "nearby period 2n" algorithm, which I posted in the first answer was my own. Another good resource is John Milnor's "Dynamic's in one complex variable"; preprint on the web, math.sunysb.edu/preprints/ims90-5.pdf . The second answer used generic complex dynamics ideas, boundary |slope|=1, plus some numerical techniques I developed to solve analytic tetration. I haven't seen the equations mapping the bulbs/carioids to the unit cirlce, except for the c=0 and c=-1 cases. – Sheldon L Sep 23 '13 at 2:34
• Anyway, I'm not claiming any original work here, it just seemed like a fun challenge that captured my imagination, and I knew I had the tools to find the Taylor series; so I did. I found this link on the web, for the period3 bulbs/cardioid. ams.org/journals/proc/1995-123-12/S0002-9939-1995-1301497-3/… For n>3 one would have to solve hopelessly complicated algebraic equations to find the fixed points; so numeric approximations would be all that is possible. I haven't seen the answer I gave posted before, but again, I'm not claiming original work. – Sheldon L Sep 23 '13 at 2:39

Edit: I answered the question of where the centers of the biggest bulbs are. I think the OP wants an approximation for the boundaries of the bulbs, which may require a more exact radius (size/2), and a better approximation of size; see my edit below.

All bulbs are periodic. If n is the period, than $f^{n}(0)=0$, where $f^{n}=(f^{n-1})^2+c$, and $f^{0}=0$. So the center of all of the bulbs (and cardiods), are all roots of algebraic equations. For n=2, you the two centers you mentioned: $c^2+c=0$, whose two roots are $c=0$ and $c=-1$. For all of these equations, you get the trivial root c=0, which can be trivially factored out. So, for n=3, after factoring out c=0, you get a cubic equation, $x^3+2x^2+x+1$, with one real root, and two complex roots. For n=4, you get a 8th order equation. Factoring out the trivial zeros, c=0 and c=-1, you get a sixth order equation, $c^6 + 3c^5 + 3c^4 + 3c^3 + 2c^2 + 1$, with two real roots, and two pairs of complex conjugate roots. This process can be continued ad nauseam, to find all values of c which are centers of bulbs. One way to find the approximate bulb size, is to say that if $f^{n}=0$, than $f^{2n}=0$ as well, but there is also another nearby zero for $f^{2n}=0$, which gives the approximate size of the bulb.

I implemented the algorithm above. Here are the results, for period(n)=2 through 19, for a size of 0.01. The size of the bulb is printed first, followed by the hyperbolic center for that bulb, which is a zero of the algebraic equation described above. The pari-gp code is posted below these results. There aren't that many big bulbs, and I didn't find any with a period>14! The algorithm uses the approximation that the size of a bulb is $1.6(z^n-z^{2n})$, where $z^{2n}$ is the nearest period 2n zero not equal to the $z^n$ zero. Edit $1.647(z^n-z^{2n})$ might be more accurate, based on the Feigenbaum constant, although that is a limiting bifurcation value, so maybe 1.61 would be better. I used an interesting routine, to estimate where the nearest 2n zero is with a cubic polynomial that factors into a quadratic, and then used Newton's method to get an exact location of the nearest 2n bulb. In the pari-gp code, if $f^n(c)$ is a zero, than "y=estim2nzero(c,n); y=invzero(y,n*2);" returns the nearest 2n zero.

You can run with smaller radius, as an optional second parameter in the periodnzero(n,l) routine. The periodnzero routine is exact, but I haven't run it for n>10 due to memory requirements. The randomperiodn routine matches the periodnzero routine results for n=5 through n=10, for a radius of 0.01, when used with 50 iterations. The randomperiodn routine only searches near the border of the main bug and the 2n bug, so it cannot find all of the smaller bulbs. For example, there is a period(6) mini-mandelbrot at -1.723, with a size of 0.0086, that the randomperiodn(6,50,0.005) doesn't find, that periodnzero(6,0.005) does find.

periodnzero(2)
0.497124226138932614 -1.00000000000000000
periodnzero(3)
0.0288243794158896630 -1.75487766624669276
0.185709338183664371 -0.122561166876653620 + 0.744861766619744237*I
periodnzero(4)
0.113351748952365738 -1.31070264133683288
0.0129942224028632946 -0.156520166833755062 + 1.03224710892283180*I
0.0868986952447806464 0.282271390766913880 + 0.530060617578525299*I
periodnzero(5)
0.0782244856281390343 -0.504340175446244000 + 0.562765761452981964*I
0.0471208943401874649 0.379513588015923745 + 0.334932305597497587*I
periodnzero(6)
0.0511889369361296355 -1.13800066665096451 + 0.240332401262098302*I
0.0282242861365737314 0.389006840569771235 + 0.215850650870819108*I
0.0433219915118591584 -0.113418655949436572 + 0.860569472501573055*I
periodnzero(7)
0.0426479226891188707 -0.622436295041293588 + 0.424878436475629184*I
0.0181617636688071151 0.376008681846767560 + 0.144749371321632865*I
0.0345136362674653694 0.121192786105906486 + 0.610611692210754212*I
periodnzero(8)
0.0246366004359986754 -1.38154748443206147
0.0294627655353046600 -0.999442387206567375 + 0.265387532468407285*I
0.0318189851964300712 -0.359102390112449160 + 0.617353453398826859*I
0.0208607361162561978 0.324819701465459679 + 0.563815622140333520*I
0.0123381391559926738 0.359031062836614452 + 0.100934876864297563*I
periodnzero(9)
0.0265592163668188400 -0.672333174859118535 + 0.337714901737740193*I
0.0177260936911318770 0.328902947296482145 + 0.410912590892970139*I
0.0185193079849810855 -0.0315529748129236492 + 0.790783175406400223*I
0.0203944122420907310 -0.210705526779016425 + 0.804635638170098617*I
randomperiodn(10,50)
0.0191394909327384273 -0.919278545244102030 + 0.247048121273991288*I
0.0203674647559298651 -1.21039963741797970 + 0.152874772077754334*I
0.0182000934472523787 0.0502718348233550644 + 0.630468552293935962*I
0.0195932658053389602 -0.533089681450721463 + 0.602309758248200881*I
0.0115562613503280188 0.408518185116081373 + 0.340038064137432468*I
randomperiodn(11,50)
0.0179922313096332877 -0.697838195122424838 + 0.279304134101365780*I
0.0101876737961521636 0.376030292897131701 + 0.266467605276236516*I
0.0142926479408153391 0.172311886095790515 + 0.570759959055939233*I
0.0170058824618214608 -0.293902905530376044 + 0.632861961802271055*I
randomperiodn(12,50)
0.0116593892097745116 -1.15193862345878364 + 0.269129812348473540*I
0.0134209900466310253 -0.871266416358976791 + 0.222317914834670266*I
0.0155654577554219361 -0.562959829953243594 + 0.471465232041119193*I
0.0116746885689485088 -1.34419625941643939 + 0.0549457739045995046*I
randomperiodn(13,50)
0.0129307773546223782 -0.712577592671273148 + 0.237792568209519750*I
0.0110764626149245331 0.0109444639810614271 + 0.638228893362735870*I
0.0129409256506876911 -0.407104083085098332 + 0.584852842868102271*I
randomperiodn(14,50)
0.0109194116271260426 -1.22997158539723289 + 0.110671419951984272*I
0.0109394099532368402 -1.05674050872754759 + 0.248827389019942445*I
0.0104598592537462953 -0.257053859166882025 + 0.639339268776504382*I
0.0111324718829115060 -0.644083434218572230 + 0.440431179466495531*I
randomperiodn(15,50)
randomperiodn(16,50)
randomperiodn(17,50)
randomperiodn(18,50)
randomperiodn(19,50)


This is the pari-gp program that generated the above results.

default(format,"g0.18");
z=1.0;
precis=precision(z);
plim=10^(-precis/2); /* precision limit */
periodnzero(j,l) = {
print("periodnzero("j")");
C0=croot(j);
for (n=1,length(C0),estimz2(C0[n],j,l));
}
/* return all the imag>0 roots of zn, iterating zn<=zn^2+x */
croot(m) = {
local(zn,Cr,n,i,v);
zn=x;
for (n=1,m-1, zn=zn^2+x; );
Cr=polroots(zn);
i=0;
v=vector(length(Cr));
for (n=1,length(Cr),
z=Cr[n];
if (((real(z)<>0) || (imag(z)<>0)) && (imag(z)>=0),
if (imag(z)==0, z=real(z));
i++;
v[i]=z;
);
);
return(v);
}
estimz2(c,n,l) = {
local(z,y,i);
z=0;
y=0;
if (l==0,l=0.01);
/* check to make sure that the z_n root is not repeating by a factor of n */
for (i=1,n-1,z=z^2+c;if (abs(z)<plim,y=i));
if (y==0,
y=estim2nzero(c,n);  /* estimate for 2n zero */
y=invzero(y,n*2); /* refinement for estimate */
/* 1.6(zn-z2n) is the approximate size of the bug or bulb */
if (abs(c-y)>l/1.6,
print(1.6*abs(c-y)" "c );
return(y);
);
return(0);
,
return(0);
);
}
/* use newton's method, centered at c */
/* iterating z=z^2+x+c with z only needing a0 and a1 terms */
invzero(c,n)={
local(z,i,a1);
z=1;
while ((abs(z)>plim) && (abs(z)<2),
z=0;
i=0;
while ((i<n) && (abs(z)<2), a1=2*a1*z+1;z=z^2+c;i++);
c=c-z/a1;
);
if (abs(z)>plim, print("invzero covergence error "c" "c+z));
return(c);
}
/* find nearest 2n zero, assuming that n is a zero */
estim2nzero(c,n) = {
local(z1,z2,y,i,a0,a1,a2,a3);
a0=0;a1=0;a2=0;a3=0;
/* iterating 2n times, centered at c, */
/* starting with z=0; z<=z^2+x+c, truncating to x^3 terms */
for (i=1,2*n,
a3=2*a3*a0+2*a1*a2;
a2=2*a2*a0+a1^2;
a1=2*a1*a0+1;
a0=a0^2+c;
);
/* assume zn=0, so z2n=0, so a0=0, */
/* divide by x, solve quadratic root closest to zero */
a0=a1;
a1=a2;
a2=a3;
y=a1^2-4*a0*a2;
z1=(-a1+sqrt(y))/(2*a2);
z2=(-a1-sqrt(y))/(2*a2);
if (abs(z1)>abs(z2),z1=z2);
return(c+z1);
}

shortperiodn(nt,nc,cnt,il) = {
local(z);
if (il==0,il=0.01,j);
z=0;for (i=1,nt,z=truncs(z^2-0+x+nc,cnt+1));
v=polroots(z);
for (n=1,cnt,
z=v[n]+nc;
y=abs(gfunc(z,nt));
if (y<0.04,
z=invzero(z,nt);
if (abs(imag(z))<1E-21, z=real(z));
if (((real(z)<>0) || (imag(z)<>0)) && (imag(z)>=0),
j=1;
while (abs(oj[j]-z)>1E-21 && (oj[j]<>0), j++);
if (oj[j]==0,
if (estimz2(z,nt,il)<>0, oj[j]=z);
);
);
);
)
}
randomperiodn(nt,n,il)= {
local(i,j);
oj=vector(128);
print("randomperiodn("nt","n")");
for (i=1,n,
shortperiodn(nt,cfromz(goldr^i),16,il); /* main cardiod */
shortperiodn(nt,-1+0.25*(goldr^i),16,il); /* period2 bulb */
);
}
/* truncate series */
truncs(z,n)={
local(y,i);
y=0;
if (n==0,n=8);
for (i=0,n-1,y=y+polcoeff(z,i)*x^i);
return(y);
}
gfunc(c,n)={local(i,z);z=0;for (i=1,n,z=z^2+c);z}
cfromz(z) = {z=z/2-1/2;z=-(z^2-1/4);}
lfromc(c) = {1/2-sqrt(1/4-c);}
gold=2/(sqrt(5)+1);
goldr=exp(2*Pi*I*gold);

• So, I wrote some code that implemented this algorithm, and it works. However, its not practical for n>8, because its apparently really memory intensive to find the roots of polynomials with more than 128 terms, although you can bump the memory on pari-gp to 1/2 a gig, and than get to 10 terms or so. So, I was going to post the roots, with radius>0.01, that I found, but the author may not be that interested. I'm exploring more random routines to find more roots for larger values of n=13, or 14 or arbitrarily large. Will post if there is interest. – Sheldon L Sep 13 '13 at 22:12
• Thanks for the insight. Is there a good article on this that is not too advanced? I am interested in the algorithm, can you post? Thanks. – PMay Sep 14 '13 at 2:31
• Posted. you could also upvote me, and give me credit for answering your question.... – Sheldon L Sep 14 '13 at 9:24
• He couldn't upvote you, having 1 rep. – Ruslan Sep 14 '13 at 9:42
• A classic book from 1990 that I own is, "Complex Dynamical Systems: The Mathematics Behind the Mandelbrot and Julia Sets", with a wonderful introductury section by Robert Devaney. You can see some of his papers online, at math.bu.edu/people/bob/papers.html – Sheldon L Sep 14 '13 at 12:24

The algebraic solutions one can find : John Stephenson : "Formulae for cycles in the Mandelbrot set", Physica A 177, 416-420 (1991); "Formulae for cycles in the Mandelbrot set II", Physica A 190, 104-116 (1992); "Formulae for cycles in the Mandelbrot set III", Physica A 190, 117-129 (1992)

Example image using this technique is here :

• Hey Adam; thanks for your reply. I didn't do much research on prior answers to how to parametrize the individual hyperbolic components. I started reading Jung's paper which looks really interesting; citeseerx.ist.psu.edu/viewdoc/… However, I may not know enough Galois theory to understand the paper. Anyway, I guess I did ok, to answer the op's question with an accurate algorithm for the Taylor series for n=3..14 for components with diameter>0.01 – Sheldon L Sep 23 '13 at 21:57
• Hi Sheldon, thanks for your work on writing the code. I have perused the literature, and I cannot find an explanation of how to numerically derive Taylor series approximations to the bulbs. Can you provide pseudo-code of your code your posted here? The Taylor series, is it a power series of angle theta, or is it in a power series of e^i*theta? – PMay Mar 10 '15 at 2:24
• – Adam Mar 12 '15 at 17:51