Fractals/Iterations in the complex plane/parabolic

"Most programs for computing Julia sets work well when the underlying dynamics is hyperbolic but experience an exponential slowdown in the parabolic case." ( Mark Braverman )

In other words it means that one can need days for making a good picture of parabolic Julia set with standard / naive algorithms.

There are 2 problems here :

• slow or lazy dynamics ( in the neighbourhood of a parabolic fixed point )
• some parts are very thin ( hard to find using standard plane scanning)

Key words

• parabolic chessboard
• parabolic implosion
• germ
• multiplicity
• Ecalle cylinders or Ecalle-Voronin cylinders ( by Jean Ecalle )
• Julia-Lavaurs sets
• petal
• The Leau-Fatou flower theorem
• Fatou coordinate
• The horn map
• Blaschke product
• Inou and Shishikura's near parabolic renormalization

eggbeater dynamics

Hand Egg beater

Here is real model of what happens in parabolic case

It is a dynamic plane for fc(z)=z^2 + 1/4. It is zoom around parabolic fixed point z=0.5. Orbits of some points inside Julia set are shown (white points)

Contunuus model of dynamics inside a 4 petals flower - dipol

Model of dynamics inside a 4 petals flower - quadrupole

Physical model : the behaviour of cake when one uses eggbeater.

The mathematical model : a 2D vector field with 2 centers ( second-order degenerate points )

The field is spinning about the centers, but does not appear to be diverging.

Petal

repelling petals around fixed point and its preimages

There is no unified definition of petals.

Petal of a flower can be :

• attracting
• repelling

Each petal is invariant under f. In other words it is mapped to itself by f.

Attracting petal P is a :

• domain (topological disc ) containing parabolic periodic point p in its boundary :

$\overline{P} \ni \left\{ p \right\}$

• trap which captures any orbit tending to parabolic point
• set contained inside component of filled-in Julia set

Petals $P_j$ are symmetric with respect to the d-1 directions :

$arg(z) = \frac{2 \Pi l}{d - 1}$

where :

• d is (to do)
• l is from 0 to d-2

Petals can have different size.

If $\lambda^n = 1$ then Julia set should approach parabolic periodic point in n directions, between n petals.

0/1

How the target set is changing along an internal ray 0

Cpp code by Wolf Jung see function parabolic from file mndlbrot.cpp ( program mandel )

To see effect :

• run Mandel
• (on parameter plane ) find parabolic point for angle 0, which is c=0.25. To do it use key c, in window input 0 and return.

C code :

  // in function uint mndlbrot::esctime(double x, double y)
  if (b == 0.0 && !drawmode && sign < 0
      && (a == 0.25 || a == -0.75)) return parabolic(x, y);
 // uint mndlbrot::parabolic(double x, double y)
 if (Zx>=0 && Zx <= 0.5 && (Zy > 0 ? Zy : -Zy)<= 0.5 - Zx) 
            { if (Zy>0) data[i]=200; // show petal
                     else data[i]=150;}

Gnuplot code :

reset
f(x,y)=  x>=0 && x<=0.5 &&  (y > 0 ? y : -y) <= 0.5 - x
unset colorbox
set isosample 300, 300
set xlabel 'x'
set ylabel 'y'
set sample 300
set pm3d map
splot [-2:2] [-2:2] f(x,y)

1/2

Cpp code by Wolf Jung see function parabolic from file mndlbrot.cpp ( program mandel ) To see effect :

• run Mandel
• (on parameter plane ) find parabolic point for angle 1/2, which is c=-0.75. To do it use key c, in window input 0 and return.

C code :

  // in function uint mndlbrot::esctime(double x, double y)
  if (b == 0.0 && !drawmode && sign < 0
      && (a == 0.25 || a == -0.75)) return parabolic(x, y);
 // uint mndlbrot::parabolic(double x, double y)
  if (A < 0 && x >= -0.5 && x <= 0 && (y > 0 ? y : -y) <= 0.3 + 0.6*x)
      {  if (j & 1) return (y > 0 ? 65282u : 65290u);
         else return (y > 0 ? 65281u : 65289u);
      }

Number of petals

In parabolic point child period coincides with parent period

For quadratic polynomials :

Multiplicity = ParentPeriod + ChildPeriod

NumberOfPetals = multiplicity - ParentPeriod

It is because in parabolic case fixed point coincidence with periodic cycle. Length of cycle ( child period) is equal to number of petals

For other polynomial maps :

• 

  One critical point and 5 petals

• 

  13 critical points and 26 petals.

• 

  14 critical points and 14 petals.

• 

  z^4-z has 3 critical points and 6 petals

f(z)	number of petals	explanation
$z^d+z$	d-1	for $z^d+z=z$ point z=0 has multiplicity d
$z^d-z$	d+2	(?)for $f^2(z)$ a root z=0 has multiplicity d+3

For f(z)= -z+z^(p+1) parabolic flower has :

• 2p petals for p odd
• p petals for p even

... ( to do )

Flower

Lea-Fatu flower

Flower with four petals and parabolic point in the center

Critical orbit for internal angle 1/5 showing 5 attracting directions

Sum of all petals creates a flower with center at parabolic periodic point.

"... an attracting petal is a set of points in a sufficient small disk around the periodic point whose forward orbits always remain in the disk under powers of return map. " ( W P Thurston : On the geometry and dynamics of Iterated rational maps)

Parabolic chessboard

Parabolic chessboard = parabolic checkerboard

See :

• Tiles: Tessellation of the Interior of Filled Julia Sets by T Kawahira
• coloured califlower by A Cheritat

Color points according to :

• the integer part of Fatou coordinate
• the sign of imaginary part

Cauliflower

Cauliflower or broccoli :

• empty ( its interior is empty ) for c outside Mandelbrot set. Julia set is a totally disconnected (
• filled cauliflower for c=1/4 on boundary of the Mandelbrot set. Julia set is a Jordan curve ( quasi circle).

• 

  c = 0 < 1/4

• 

  c = 0.25 = 1/4 t he root of the main cardioid. Julia set is a cauliflower

• 

  c = 0.255 > 1/4 "imploded cauliflower"

• 

  c = 0.258 > 1/4

• 

  c = [0.285, 0.01]

Pleae note that :

• size of image differs because of different z-planes.
• different algorithms are used so colours are hard to compare

Bifurcation of the Cauliflower

How Julia set changes along real axis ( going from c=0 thru c=1/4 and futher ) :

Perturbation of a function $f(z)$ by complex $\epsilon$ :

$g(z) = f(z) + \epsilon$

When one add epsilon > 0 ( move along real axis toward + infinity ) there is a bifurcation of parabolic fixed point :

• attracting fixed point ( epsilon<0 )
• one parabolic fixed point ( epsilon = 0 )
• one parabolic fixed point splits up into two conjugate repelling fixed points ( epsilon > 0 )

"If we slightly perturb with epsilon<0 then the parabolic fixed point splits up into two real fixed points on the real axis (one attracting, one repelling). "

See :

• demo 2 page 9 in program Mandel by Wolf Jung

parabolic implosion

Video on YouTube

How to compute parabolic c values

Description

Period 1

One can easily compute boundary point c

$c = c_x + c_y*i$

of period 1 hyperbolic component ( main cardioid) for given internal angle ( rotation number) t using this cpp code by Wolf Jung

t *= (2*PI); // from turns to radians
cx = 0.5*cos(t) - 0.25*cos(2*t); 
cy = 0.5*sin(t) - 0.25*sin(2*t); 

or this Maxima CAS code :

 
/* conformal map  from circle to cardioid ( boundary
 of period 1 component of Mandelbrot set */
F(w):=w/2-w*w/4;

/* 
circle D={w:abs(w)=1 } where w=l(t,r) 
t is angle in turns ; 1 turn = 360 degree = 2*Pi radians 
r is a radius 
*/
ToCircle(t,r):=r*%e^(%i*t*2*%pi);

GiveC(angle,radius):=
(
 [w],
 /* point of  unit circle   w:l(internalAngle,internalRadius); */
 w:ToCircle(angle,radius),  /* point of circle */
 float(rectform(F(w)))    /* point on boundary of period 1 component of Mandelbrot set */
)$

compile(all)$

/* ---------- global constants & var ---------------------------*/
Numerator :1;
DenominatorMax :10;
InternalRadius:1;

/* --------- main -------------- */
for Denominator:1 thru DenominatorMax step 1 do
(
 InternalAngle: Numerator/Denominator,
 c: GiveC(InternalAngle,InternalRadius),
 display(Denominator),
 display(c),
  /* compute fixed point */
 alfa:float(rectform((1-sqrt(1-4*c))/2)), /* alfa fixed point */
 display(alfa)
 )$

Internal angle t	c	fixed point $z_{alfa}$
1/1	0.25	0.5
1/2	-0.75	-0.5
1/3	0.64951905283833*%i-0.125	0.43301270189222*%i-0.25
1/4	0.5*%i+0.25	0.5*%i
1/5	0.32858194507446*%i+0.35676274578121	0.47552825814758*%i+0.15450849718747
1/6	0.21650635094611*%i+0.375	0.43301270189222*%i+0.25
1/7	0.14718376318856*%i+0.36737513441845	0.39091574123401*%i+0.31174490092937
1/8	0.10355339059327*%i+0.35355339059327	0.35355339059327*%i+0.35355339059327
1/9	0.075191866590218*%i+0.33961017714276	0.32139380484327*%i+0.38302222155949
1/10	0.056128497072448*%i+0.32725424859374	0.29389262614624*%i+0.40450849718747

For internal angle n/p parabolic period p cycle consist of one z-point with multiplicity p and multiplier = 1.0 . This point z is equal to fixed point $z_{alfa}$

Period 2

// cpp code by W Jung http://www.mndynamics.com

t *= (2*PI);  // from turns to radians
cx = 0.25*cos(t) - 1.0;
cy = 0.25*sin(t);  

Periods 1-6

/* 


batch file for Maxima CAS 
computing bifurcation points for period 1-6

 Formulae for cycles in the Mandelbrot set II
Stephenson, John; Ridgway, Douglas T.
Physica A, Volume 190, Issue 1-2, p. 104-116.
*/

kill(all);
remvalue(all);


start:elapsed_run_time ();

/* ------------ functions ----------------------*/

/* exponential for of complex number with angle in turns */
 /* "exponential form prevents allroots from working", code by Robert P. Munafo */ 

GivePoint(Radius,t):=rectform(ev(Radius*%e^(%i*t*2*%pi), numer))$ /* gives point of unit circle for angle t in turns */

GiveCirclePoint(t):=rectform(ev(%e^(%i*t*2*%pi), numer))$ /* gives point of unit circle for angle t in turns Radius = 1 */

/* gives a list of iMax points of unit circle */
GiveCirclePoints(iMax):=block(
 [circle_angles,CirclePoints],
 CirclePoints:[],
 circle_angles:makelist(i/iMax,i,0,iMax),
 for t in circle_angles do CirclePoints:cons(GivePoint(1,t),CirclePoints),
 return(CirclePoints) /* multipliers */
)$

/* http://commons.wikimedia.org/wiki/File:Mandelbrot_set_Components.jpg 
Boundary equation  b_n(c,P)=0 
    defines relations between hyperbolic components and unit circle for given period n ,
    allows computation of exact coordinates of hyperbolic componenets.

b_n(w,c), is boundary polynomial ( implicit function of 2 variables ).

*/

GiveBoundaryEq(P,n):=
block(
 if n=1 then return(c + P^2 - P),
 if n=2 then return(- c + P - 1),
 if n=3 then return(c^3 + 2*c^2 - (P-1)*c + (P-1)^2),
 if n=4 then return( c^6 + 3*c^5 + (P+3)* c^4 + (P+3)* c^3  - (P+2)*(P-1)*c^2 - (P-1)^3),
 if n=5 then return(c^15 + 8*c^14 + 28*c^13 + (P + 60)*c^12 + (7*P + 94)*c^11 + 
  (3*P^2 + 20*P + 116)*c^10 + (11*P^2 + 33*P + 114)*c^9 + (6*P^2 + 40*P + 94)*c^8 + 
  (2*P^3 - 20*P^2 + 37*P + 69)*c^7 + (3*P - 11)*(3*P^2 - 3*P - 4)*c^6 + (P - 1)*(3*P^3 + 20*P^2 - 33*P - 26)*c^5 +
  (3*P^2 + 27*P + 14)*(P - 1)^2*c^4 - (6*P + 5)*(P - 1)^3*c^3 + (P + 2)*(P - 1)^4*c^2 - c*(P - 1)^5  + (P - 1)^6),
if n=6 then return( c^27+
13*c^26+
78*c^25+
(293 - P)*c^24+
(792 - 10*P)*c^23+
(1672 - 41*P)*c^22+
(2892 - 84*P - 4*P^2)*c^21+
(4219 - 60*P - 30*P^2)*c^20+
(5313 + 155*P - 80*P^2)*c^19+
(5892 + 642*P - 57*P^2 + 4*P^3)*c^18+
(5843 + 1347*P + 195*P^2 + 22*P^3)*c^17+
(5258 + 2036*P + 734*P^2 + 22*P^3)*c^16+
(4346 + 2455*P + 1441*P^2 - 112*P^3 + 6*P^4)*c^15 + 
(3310 + 2522*P + 1941*P^2 - 441*P^3 + 20*P^4)*c^14 + 
(2331 + 2272*P + 1881*P^2 - 853*P^3 - 15*P^4)*c^13 + 
(1525 + 1842*P + 1344*P^2 - 1157*P^3 - 124*P^4 - 6*P^5)*c^12 + 
(927 + 1385*P + 570*P^2 - 1143*P^3 - 189*P^4 - 14*P^5)*c^11 + 
(536 + 923*P - 126*P^2 - 774*P^3 - 186*P^4 + 11*P^5)*c^10 + 
(298 + 834*P + 367*P^2 + 45*P^3 - 4*P^4 + 4*P^5)*(1-P)*c^9 + 
(155 + 445*P - 148*P^2 - 109*P^3 + 103*P^4 + 2*P^5)*(1-P)*c^8 + 
2*(38 + 142*P - 37*P^2 - 62*P^3 + 17*P^4)*(1-P)^2*c^7 + 
(35 + 166*P + 18*P^2 - 75*P^3 - 4*P^4)*((1-P)^3)*c^6 + 
(17 + 94*P + 62*P^2 + 2*P^3)*((1-P)^4)*c^5 + 
(7 + 34*P + 8*P^2)*((1-P)^5)*c^4 + 
(3 + 10*P + P^2)*((1-P)^6)*c^3 + 
(1 + P)*((1-P)^7)*c^2 +
-c*((1-P)^8) + (1-P)^9)
)$


/* gives a list of points c on boundaries on all components for give period */
GiveBoundaryPoints(period,Circle_Points):=block(
 [Boundary,P,eq,roots],
  Boundary:[],
 for m in Circle_Points do (/* map from reference plane to parameter plane */
  P:m/2^period,
  eq:GiveBoundaryEq(P,period), /* Boundary equation  b_n(c,P)=0  */
  roots:bfallroots(%i*eq),
  roots:map(rhs,roots),
  for root in roots do Boundary:cons(root,Boundary)),
  return(Boundary)
)$


/* divide llist of roots to 3 sublists = 3  components */
/* ---- boundaries of period 3 components 
period:3$
Boundary3Left:[]$
Boundary3Up:[]$
Boundary3Down:[]$

Radius:1;

 for m in CirclePoints do (
  P:m/2^period,
  eq:GiveBoundaryEq(P,period),
  roots:bfallroots(%i*eq),
  roots:map(rhs,roots),
  for root in roots do 
     (
       if realpart(root)<-1  then Boundary3Left:cons(root,Boundary3Left),
       if (realpart(root)>-1 and imagpart(root)>0.5) 
            then Boundary3Up:cons(root,Boundary3Up),
       if (realpart(root)>-1 and imagpart(root)<0.5) 
            then Boundary3Down:cons(root,Boundary3Down)
               
     )

)$
--------- */


/* gives a list of parabolic points for given : period and internal angle */
GiveParabolicPoints(period,t):=block
(
 [m,ParabolicPoints,P,eq,roots],
 m: GiveCirclePoint(t), /* root of unit circle, Radius=1, angle t=0 */
 ParabolicPoints:[],
 /* map from reference plane to parameter plane */
 P:m/2^period,
 eq:GiveBoundaryEq(P,period), /* Boundary equation  b_n(c,P)=0  */
 roots:bfallroots(%i*eq),
 roots:map(rhs,roots),
 for root in roots do ParabolicPoints:cons(float(root),ParabolicPoints),
 return(ParabolicPoints) 

)$


compile(all)$



/* ------------- constant values ----------------------*/

fpprec:16; 





/* ------------unit circle on a w-plane -----------------------------------------*/
a:GiveParabolicPoints(6,1/3);
a$

How to draw parabolic Julia set

• Three methods of drawing
• 

  Checking iteration

• 

  modified DEM

• 

  checking angle

All points of interior of filled Julia set tend to one periodic orbit ( or fixed point ). This point is in Julia set and is weakly attracting. One can analyse only behevior near parabolic fixed point. It can be done using critical orbits.

There are two cases here : easy and hard.

If the Julia set near parabolic fixed point is like n-th arm star ( not twisted) then one can simply check argument of of zn, relative to the fixed point. See for example z+z^5. This is an easy case.

In the hard case Julia set is twisted around fixed.

Estimation from exterior

Escape time

Description

Long iteration method

Long iteration method

Dynamic rays

One can use periodic dynamic rays landing on parabolic fixed point to find narrow parts of exterior.

Let's check how many backward iterations needs point on periodic ray with external radius = 4 to reach distance 0.003 from parabolic fixed point :

period	Inverse iterations	time
1	340	0m0.021s
2	55 573	0m5.517s
3	8 084 815	13m13.800s
4	1 059 839 105	1724m28.990s

Estimation from interior

Julia set is a boundary of filled-in Julia set Kc.

• find points of interior of Kc
• find boundary of interior of Kc using edge detection

If components of interior are lying very close to each other then find components using :

color = LastIteration % period

For parabolic components between parent and child component :

periodOfChild = denominator*periodOfParent  
color = iLastIteration % periodOfChild 

where denominator is a denominator of internal angle of parent comonent of Mandelbrot set.

Angle

"if the iterate zn of tends to a fixed parabolic point, then the initial seed z0 is classified according to the argument of zn−z0, the classification being provided by the flower theorem " ( Mark McClure )

Attraction time

Various types of dynamics

Attraction time for various denominators

Interior of filled Julia set consist of components. All comonents are preperiodic, some of them are periodic ( immediate basin of attraction).

In other words :

• one iteration moves z to another component ( and whole component to another component)
• all point of components have the same attraction time ( number of iteration needed to reach target set around attractor)

It is possible to use it to color components. Because in parabolic case attractor is weak ( weakly attracting) it needs a lot of iterations for some points to reach it.

Here are some example values :

iWidth  = 1001 // width of image in pixels
PixelWidth  = 0.003996  
AR  = 0.003996 // Radius around attractor
denominator  = 1 ; Cx  = 0.250000000000000; Cy  = 0.000000000000000 ax  = 0.500000000000000; ay  = 0.000000000000000   
denominator  = 2 ; Cx  = -0.750000000000000; Cy  = 0.000000000000000 ax  = -0.500000000000000; ay  = 0.000000000000000   
denominator  = 3 ; Cx  = -0.125000000000000; Cy  = 0.649519052838329 ax  = -0.250000000000000; ay  = 0.433012701892219  
denominator  = 4 ; Cx  = 0.250000000000000; Cy  = 0.500000000000000 ax  = 0.000000000000000; ay  = 0.500000000000000   
denominator  = 5 ; Cx  = 0.356762745781211; Cy  = 0.328581945074458 ax  = 0.154508497187474; ay  = 0.475528258147577   
denominator  = 6 ; Cx  = 0.375000000000000; Cy  = 0.216506350946110 ax  = 0.250000000000000; ay  = 0.433012701892219    

denominator  = 1 ;   i =               243.000000 
denominator  = 2 ;   i =            31 171.000000 
denominator  = 3 ;   i =         3 400 099.000000 
denominator  = 4 ;   i =       333 293 206.000000 
denominator  = 5 ;   i =    29 519 565 177.000000 
denominator  = 6 ;   i = 2 384 557 783 634.000000 

where :

C = Cx + Cy*i 
a = ax + ay*i // fixed point alpha
i // number of iterations after which critical point z=0.0 reaches disc around fixed point alpha with radius AR
denominator of internal angle ( in turns )
internal angle =  1/denominator

Note that attraction time i is proportional to denominator.

Now you see what means weakly attracting.

One can :

• use brutal force method ( Attracting radius < pixelSize; iteration Max big enough to let all points from interior reach target set; long time or fast computer )
• find better method (:-)) if time is to long for you

Interior distance estimation

Trap

Estimation from interior and exterior

Julia set is a common boundary of filled-in Julia set and basin of attraction of infinity.

• find points of interior/components of Kc
• find escaping points
• find boundary points using Sobel filter

It works for denominator up to 4.

Inverse iteration of repelling points

Inverse iteration of alfa fixed point. It works good only for cuting point ( where external rays land). Other points still are not hitten.

Bof61

• Using bof61 for coloring interior
• example

Gallery

• Critical orbits in case of one critical point
• 

  critical orbits for internal angle from 1/1 to 1/10. True attracting directions

• 

  N-th arm stars for n from 1 to 10. Schematic attracting directions

• 

  critical orbit, attracting and repelling vectors for internal angle 1/3

• from period 1 thru ...
• 

  thru internal ray 1/1; c = 0.25 = 1/4 the root of the main cardioid. Julia set is a cauliflower

• 

  from period 1 to 2 thru internal ray 1/2 ; c=-0.75 Julia set is a San Marco fractal

• 

  from period 1 thru internal ray 1/3 ; Juia set is a Douady fat rabbit

• 

  internal angle 1/4

• 

  internal angle 1/5

• 

  internal angle 1/7

• from period 2 thru ...
• 

  from 2 thru internal ray 1/1 ; c=-0.75 Julia set is a San Marco fractal

• 

  from period 2 thru 1/2

• 

  from period 2 thru 1/3

• 

  from period 2 thru 1/4

• from period 3 thru ...
• 

  from period 3 thru 1/3

For other polynomial maps see here

See also

• Image : Nonstandard Parabolic by Cheritat
• Julia set of parabolic case in Maxima CAS
• The parabolic Mandelbrot set Pascale ROESCH (joint work with C. L. PETERSEN
• PARABOLIC IMPLOSION A MINI-COURSE by ArnaudCheritat
• Workshop on parabolic implosion 2010
• The renormalization for parabolic fixed points... Mitsushiro Shishikura and Hiroyuki Inou

References