Combined Methods of Depicting Julia Sets

 

Fig: Quadratic Julia sets depicted by combined methods. (a) A Siegel disc (). Exterior (blue) by level set, Julia (white) by modified inverse iteration, irrational flow (green) by the sine of the velocity. (b) Super-attracting period 3 (Solution of ). Exterior by continuous potential, Julia by modified inverse iteration, interior basins by level set. (c) Seahorse valley dendrite () by distance estimator, exterior by level set. (d) Parabolic bifurcation period1 to 20 ( adjacent to ). Exterior by level set, Julia by modified inverse iteration, petals by velocity.

Level Set Method LSM

Begin with the point for a fixed c and iterate . If we are depicting the Mandelbrot parameter plane, we begin instead with the critical point and iterate for each . If the point escapes a circle, say of radius 10, , we colour it by the number of iterations.  If it remains bounded after a fixed number of iterations, we assume it cannot escape and colour it black (or white in the above examples).

 

This method will work for the Julia set of any function provided we can determine basins of attraction of fixed or periodic points, to apply the algorithm to. In the case of  this is  and the Julia set is the boundary of the basin of attraction of . In other functions it may be finite and multiple. In the case of transcendental functions, such as cos and exp,  is an essential singularity and the Julia set is the closure of its basin of attraction and the target set may be  or . This method will also work for the parameter plane (Mandelbrot set) of any function provided we can locate and handle the critical points and establish individual parameter planes for each.

 

Variation 1 Internal Basins: 

If the Julia set has an internal attracting periodic point we can test for this by also finding the periodic cycle by first iterating the critical point until it becomes within  of being periodic and then colour a non-escaping point by the number of iterations to bring  within  of a point on the cycle.

 

Variation 2 Continuous Potential (Quadratics only):

We can derive a continuous potential with an even gradation but less sensitivity at the boundary of  by using the continuous potential formula  for an escaping point.

Variation 3 Distance Estimator (Quadratics only):

For dealing with dendritic Julia sets, or the Mandelbrot set, to demonstrate connectivity, we can make two additional tests to detect proximity to :

(a) Test for overflow of the derivative of any point on the escaping orbit:

(b) Test for distance

We can calculate orbit derivatives iteratively as follows  for the Julia set and  for the Mandelbrot set, since we are differentiating with respect to c rather than z. In practice, O ~ 1600000,  ~ 0.1. The derivative overflow tests for highly repulsive dynamics adjacent to .

 

Variation 4 Discrete Velocity of non-attracting Basins and Petals:

Compute, for the points that don?t escape, the average discrete velocity  on the orbit.

 

Herman ring  Irrational annular flow depicted

 by sine of the discrete velocity, interior basin level sets red, exterior blue.

Inverse Iteration Method IIM

First find a repelling fixed point , , by solving . For , one of the two fixed points is always a repeller as , unless . Now plot the two inverse images  of this point and repeat to form the   n-th inverse iterates. This method requires a heap, or some equivalent data structure, to keep track of the branching tree of inverse iterates. If memory is exceeded we can randomly plot one or other roots and its pre-images. However this method has the problem that it is computationally intractable because the points are exponentially unevenly distributed over the Julia set, due to multifractality (fractal redistribution of the probabilities), resulting in the inverse mapping being strongly contractive to some features leaving others unrepresented, and thus fails to represent significant features of , even with exhaustive computation times.

 

Variation Modified Inverse Iteration MIIM:

We cut off the sub tree from a given  if the derivative .

This eliminates dominant highly contractive regions of the inverse iteration, which have already been registered. We can calculate successive derivatives iteratively .

Sample C code for heap processing of MIIM using pointers

 


typedef double cno[2];

 

const short finish = 10000;

const short scal = 100;

const double gam = 2.0;

const double huge = 10000;

 

struct fork

{

fork*        back;

fork*        side;

fork*        fwd;

cno           num;

double     drv;

};

 

double a,b;

int      nw, level,mysiz;

cno  root1, root2;

fork* here;

Handle    mSet, mMat;

long iters;

Rect                                             boundsRect;

WindowPtr                                 theWindow;

short globalRef,globaliRef;

long mapSize, PackScreen(Ptr);

 

cno* csq(cno cc)   //gives 1 of 2 complex square roots - the other is symmetric

{

double   rt;

cno csqrt;

double sq2=sqrt(2);

 

  if (cc[1] != 0)

  {

    rt = sqrt(cc[0] + sqrt(cc[0] * cc[0] + cc[1] * cc[1]));

    csqrt[0] = rt / sq2;

                  csqrt[1] = cc[1] / sq2 / rt;

   }

  else if (cc[0] >= 0)

   {

    csqrt[0] = sqrt(cc[0]);

    csqrt[1] = 0;

   }

  else

   {

    csqrt[1] = sqrt(-cc[0]);

    csqrt[0] = 0;

   }

return(&csqrt);

}

 

void iterate(cno val)

{

short i;

cno* myp;

val[0]=val[0]-a;

val[1]=val[1]-b;

myp=csq(val);

for(i=0;i<2;i++) root1[i]=(*myp)[i];

for(i=0;i<2;i++) root2[i]=-(*myp)[i];

}

 

int okay(void)

{

if(here->drv<huge) return(1);

 else return(0);

}

 

void fwdispold (void)

{

nw=0;

level+=-1;

here=here->fwd;

free(here->back);

}

 

void sideptr(cno* myp)

{

short i;

here->side=(fork*)calloc(1,sizeof(fork));

here->side->fwd=here;

for(i=0;i<2;i++) here->side->num[i]=(*myp)[i];

here->side->drv=gam*here->drv*sqrt((*myp)[0]*(*myp)[0]+(*myp)[1]*(*myp)[1]);

}

 

void backptr(cno* myp)

{

short i;

nw=1;

here->back=(fork*)calloc(1,sizeof(fork));

here->back->fwd=here;

for(i=0;i<2;i++) here->back->num[i]=(*myp)[i];

here->back->drv=gam*here->drv*sqrt((*myp)[0]*(*myp)[0]+(*myp)[1]*(*myp)[1]);

here=here->back;

level+=1;

}

 

void plot(cno* myp)

{

cno val;

short i;

for(i=0;i<2;i++) val[i]=scal*(1.5+(*myp)[i]);

MoveTo(val[0],val[1]);

LineTo(val[0],val[1]);

}

 

void perform(void)

{

short i;

cno  c, z, *zz;

 

a=0.27334;

b=0.00742;

 c[0] = a;

 c[1] = b;

 z[0] = 1 - 4 * c[0];

 z[1] = -4 * c[1];

 zz = csq(z);

 z[0] = (1 - (*zz)[0]) / 2;

 z[1] = -(*zz)[1] / 2;

 if (z[0] * z[0] + z[1] * z[1] < 0.5)

  {

   z[0] = (1 + (*zz)[0]) / 2;

   z[1] = (*zz)[1] / 2;    //find the greatest i.e. repelling root

  }

nw=1;

level=0;

here=(fork*)calloc(1,sizeof(fork));

for(i=0;i<2;i++) here->num[i]=z[i];

here->drv=gam*sqrt(z[0]*z[0]+z[1]*z[1]);

do

{

if(nw)

                  {

                  if(okay())

                                    {

                                    plot(&(here->num));

                                    if(level<finish)

                                                      {

                                                      iterate(here->num);

                                                      sideptr(&root1);

                                                      backptr(&root2);

                                                      }

                                    else fwdispold();

                                    }

                  else fwdispold();

                  }

else

                  if(here->side!=0)

                                    {

                                    //backptr(&(here->side->num));

                                    here->back=here->side;

                                    here->side=0;

                                    here=here->back;

                                    level+=1;

                                    //free(here->fwd->side);

                                    nw=1;

                                    }

                  else fwdispold();

iters=iters+1;

}

while(level>=0);

}