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.

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.

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 .

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);

}