Fractals/Iterations in the complex plane/Julia set

This book shows how to code different algorithms for drawing sets in dynamical plane : Julia, Filled-in Julia or Fatou sets for complex quadratic polynomial. It is divided in 2 parts :

• description of various algorithms
• descriptions of technics for visualisation of various sets in dynamic plane
  • Julia set
  • Fatou set
    • basin of attraction of infinity ( open set)
    • basin of attraction of finite attractor

Various types of dynamics needs various algorithms

Algorithms

Methods based on speed of attraction

Here color is proportional to speed of attraction ( convergence to attractor). These methods are used in Fatou set.

Basin of attraction to infinity = exterior of filled-in Julia set and The Divergence Scheme = Escape Time Method ( ETM )

First read definitions

Here one computes forward iterations of a complex point Z₀:

$Z_{n+1} = Z_n^2 + C$

Here is function which computes the last iteration, that is the first iteration that lands in the target set ( for example leaves a circle around the origin with a given escape radius ER ) for the iteration of the complex quadratic polynomial above. It is a iteration ( integer) for which (abs(Z)>ER).

C version ( here ER2=ER*ER) using double floating point numbers ( without complex type numbers) :

 int GiveLastIteration(double Zx, double Zy, double Cx, double Cy, int IterationMax, int ER2)
 {
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  int i=0;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  while (i<IterationMax && (Zx2+Zy2<ER2) ) /* ER2=ER*ER */
  {
   Zy=2*Zx*Zy + Cy;
   Zx=Zx2-Zy2 +Cx;
   Zx2=Zx*Zx;
   Zy2=Zy*Zy;
   i+=1;
  }
  return i;
 }

C++ versions:

 int GiveLastIteration(complex C,complex Z , int imax, int ER)
  {
   int i; // iteration number
   for(i=0;i<=imax-1;i++) // forward iteration
    {
      Z=Z*Z+C; // overloading of operators
      if(abs(Z)>ER) break;
    }
   return i;
 }

#include <complex> // C++ complex library 
 
// bailout2 = bailout * bailout
// this function is based on function esctime from mndlbrot.cpp 
// from program mandel ver. 5.3 by Wolf Jung
// http://www.mndynamics.com/indexp.html 
 
int escape_time(complex<double> Z, complex<double> C , int iter_max,  double bailout2)
{ 
  // z= x+ y*i   z0=0
  long double x =Z.real(), y =Z.imag(),  u ,  v ;
  int iter; // iteration
  for ( iter = 0; iter <= iter_max-1; iter++)
  { u = x*x; 
    v = y*y;
    if ( u + v <= bailout2 ) 
       { 
         y = 2 * x * y + C.imag();
         x = u - v + C.real(); 
       } // if
    else break;
  } // for 
  return iter;
} // escape_time

Delphi version ( using user defined complex type, cabs and f functions )

function GiveLastIteration(z,c:Complex;ER:real;iMax:integer):integer;
  var i:integer;
  begin
  i:=0;
  while (cabs(z)<ER) and (i<iMax) do
    begin
      z:= f(z,c);
      inc(i);
    end;
  result := i;
end;

where :

type complex = record x, y: real; end;
 
function cabs(z:complex):real;
begin
  cabs:=sqrt(z.x*z.x+z.y*z.y)
end;
 
function f(z,c:complex):complex; // complex quadratic polynomial
  var tmp:complex;
begin
  tmp.x := (z.x*z.x) - (z.y*z.y) + c.x;
  tmp.y := 2*z.x*z.y + c.y ;
  result := tmp;
 
end;

Delphi version without explicit definition of complex numbers :

function GiveLastIteration(zx0,zy0,cx,cy,ER2:extended;iMax:integer):integer;
  // iteration of z=zx+zy*i under fc(z)=z*z+c
  // where c=cx+cy*i
  // until abs(z)<ER  ( ER2=ER*ER )  or i>=iMax
  var i:integer;
      zx,zy,
      zx2,zy2:extended;
  begin
  zx:=zx0;
  zy:=zy0;
  zx2:=zx*zx;
  zy2:=zy*zy;
 
  i:=0;
  while (zx2+zy2<ER2) and (i<iMax) do
    begin
      zy:=2*zx*zy + cy;
      zx:=zx2-zy2 +cx;
      zx2:=zx*zx;
      zy2:=zy*zy;
      //
      inc(i);
    end;
  result := i;
end;

Euler version by R. Grothmann ( with small change : from z^2-c to z^2+c)

function iter (z,c,n=100) ...

h=z;
loop 1 to n;
h=h^2 + c;
if totalmax(abs(h))>1e20; m=#; break; endif;
end;
return {h,m};
endfunction

Lisp version

This version uses complex numbers. It makes the code short but is also inefficien.

((DEFUN GIVELASTITERATION (Z_0 _C IMAX ESCAPE_RADIUS) 
  (SETQ Z Z_0) 
  (SETQ I 0)
  (LOOP WHILE (AND (< I IMAX) (< (ABS Z) ESCAPE_RADIUS)) DO 
    (INCF I)
    (SETQ Z (+ (* Z Z) _C)))
   I)

Maxima version :

/* easy to read but very slow version, uses complex type numbers */ 
GiveLastIteration(z,c):=
 block([i:0],
  while abs(z)<ER and i<iMax
     do (z:z*z + c,i:i+1),
  i)$

/* faster version, without use of complex type numbers,
   compare with c version, ER2=ER*ER */  
GiveLastIter(zx,zy,cx,cy,ER2,iMax):=
block(
 [i:0,zx2,zy2],
 zx2:zx*zx,
 zy2:zy*zy,
 while (zx2+zy2<ER2) and i<iMax do
 (      
  zy:2*zx*zy + cy,
  zx:zx2-zy2 +cx,
  zx2:zx*zx,
  zy2:zy*zy,
  i:i+1
 ),
 return(i)
);

Boolean Escape time

Algorithm: for every point z of dynamical plane (z-plane) compute iteration number ( last iteration) for which magnitude of z is greater than escape radius. If last_iteration=max_iteration then point is in filled-in Julia set, else it is in its complement (attractive basin of infinity ). Here one has 2 options, so it is named boolean algorithm.

if (LastIteration==IterationMax)
   then color=BLACK;   /* bounded orbits = Filled-in Julia set */
   else color=WHITE;  /* unbounded orbits = exterior of Filled-in Julia set  */

In theory this method is for drawing Filled-in Julia set and its complement ( exterior), but when c is Misiurewicz point ( Filled-in Julia set has no interior) this method draws nothing. For example for c=i . It means that it is good for drawing interior of Filled-in Julia set .

ASCII graphic

; common lisp
(loop for y from -2 to 2 by 0.05 do
      (loop for x from -2 to 2 by 0.025 do
                (let* ((z (complex x y))
                        (c (complex -1 0))
                        (iMax 20)
                        (i 0))
 
                (loop      while (< i iMax ) do 
                        (setq z (+ (* z z) c))
                        (incf i)
                        (when (> (abs z) 2) (return i)))
 
 
 
           (if (= i iMax) (princ (code-char 42)) (princ (code-char 32)))))
      (format t "~%"))

PPM file with raster graphic

Filled-in Julia set for c= =-1+0.1*i. Image and C source code

Integer escape time = Level Sets of the Basin of Attraction of Infinity = Level Sets Method= LSM/J

level sets of escape time for C=0 and 0.9<Z<1.5

Escape time measures time of escaping to infinity ( infinity is superattracting point for polynomials). Time is measured in steps ( iterations = i) needed to escape from circle of given radius ( ER= Escape Radius).

One can see few things:

• this is discontinuous function
• i is iMax for z in Filled-in Julia set
• i=0 for x0>ER
• this is nonlinear function

Level sets here are sets of points with the same escape time. Here is algorithm of choosing color in black & white version.

  if (LastIteration==IterationMax)
   then color=BLACK;   /* bounded orbits = Filled-in Julia set */
   else   /* unbounded orbits = exterior of Filled-in Julia set  */
        if ((LastIteration%2)==0) /* odd number */
           then color=BLACK; 
           else color=WHITE;

Normalized iteration count (real escape time or fractional iteration)

Escape time for C=0 and 0.5<Z<2.5

Math formula :

$\nu(z)=\lim\limits_{i\to\infty}( i-\log_2\log_2|z_i|)\ .$

Maxima version :

GiveNormalizedIteration(z,c,E_R,i_Max):=
/* */
block(
 [i:0,r],
 while abs(z)<E_R and i<i_Max
   do (z:z*z + c,i:i+1),
 r:i-log2(log2(cabs(z))),
 return(float(r))
)$

In Maxima log(x) is a natural (base e) logarithm of x. To compute log2 use :

log2(x) := log(x) / log(2);

Level Curves of escape time Method = eLCM/J

• 

  edge detection Image and C code

• 

  preimages of circle

These curves are boundaries of Level Sets of escape time ( eLSM/J ). They can be drawn using these methods:

• edge detection of Level Curves ( =boundaries of Level sets).
  • Algorithm based on paper by M. Romera et al.
  • Sobel filter
• drawing lemniscates = curves $L_n=\{z: abs(z_n)=ER \}\,$, see explanation and source code
• drawing circle $L_0=\{z: abs(z)=ER \}\,$ and its preimages. See this image, explanation and source code
• method described by Harold V. McIntosh

/* Maxima code : draws lemniscates of Julia set */
c: 1*%i;
ER:2;
z:x+y*%i;
f[n](z) := if n=0 then z else (f[n-1](z)^2 + c);
load(implicit_plot); /* package by Andrej Vodopivec */ 
ip_grid:[100,100];
ip_grid_in:[15,15];
implicit_plot(makelist(abs(ev(f[n](z)))=ER,n,1,4),[x,-2.5,2.5],[y,-2.5,2.5]);

Basin of attraction of finite attractor = interior of filled-in Julia set

• How to find periodic attractor ?
• How many iterations is needed to reach attractor ?

Distance between points and iteration

Components of Interior of Filled Julia set ( Fatou set)

• use limited color ( palette = list of numbered colors)
• find period of attracting cycle
• find one point of attracting cycle
• compute number of iteration after when point reaches the attractor
• color of component=iteration % period
• use edge detection for drawing Julia set

• 

  Components of Fatou set for period 4

• 

  Components of Fatou set for period 3

• 

  In case of Siegel disc critical orbit is a boundary Siegel disc compponent. All other componnats are preimages of this component

Internal Level Sets

See :

• algorithm 0 of program Mandel by Wolf Jung

• 

  c = 0.5*i. Image and source code

• 

  c = -0.584895. Image and source code

Decomposition of target set

• Decomposition
• 

  Binary decomposition of whole dynamical plane with circle Julia set c = 0

• 

  Modified binary decomposition of whole dynamical plane with circle Julia set. Image with src code

• 

  Another modification of binary decomposition

• 

  Modified decomposition of dynamical plane with dendrite Julia set. Image with src code

• 

  Modified decomposition of dynamical plane with basilica Julia set. Image with src code

Binary decomposition

Here color of pixel ( exterior of Julia set) is proportional to sign of imaginary part of last iteration .

Main loop is the same as in escape time.

In other words target set is decompositioned in 2 parts ( binary decomposition) :

$T_b+=\{z: abs(z_n) > ER ~~\mbox{and}~~ im(z_n)>0 \} \,$

$T_b-=\{z: abs(z_n) > ER ~~\mbox{and}~~ im(z_n)<=0 \} \,$

Algorithm in pseudocode ( Im(Zn) = Zy ) :

if (LastIteration==IterationMax)
   then color=BLACK;   /* bounded orbits = Filled-in Julia set */
   else   /* unbounded orbits = exterior of Filled-in Julia set  */
        if (Zy>0) /* Zy=Im(Z) */
           then color=BLACK; 
           else color=WHITE;

Modified decomposition

Here exterior of Julia set is decompositioned into radial level sets.

It is because main loop is without bailout test and number of iterations ( iteration max) is constant.

It creates radial level sets.

  for (Iteration=0;Iteration<8;Iteration++)
 /* modified loop without checking of abs(zn) and low iteration max */
    {
    Zy=2*Zx*Zy + Cy;
    Zx=Zx2-Zy2 +Cx;
    Zx2=Zx*Zx;
    Zy2=Zy*Zy;
   };
  iTemp=((iYmax-iY-1)*iXmax+iX)*3;        
  /* --------------- compute  pixel color (24 bit = 3 bajts) */
 /* exterior of Filled-in Julia set  */
 /* binary decomposition  */
  if (Zy>0 ) 
   { 
    array[iTemp]=255; /* Red*/
    array[iTemp+1]=255;  /* Green */ 
    array[iTemp+2]=255;/* Blue */
  }
  if (Zy<0 )
   {
    array[iTemp]=0; /* Red*/
    array[iTemp+1]=0;  /* Green */ 
    array[iTemp+2]=0;/* Blue */    
  };

It is also related with automorphic function for the group of Mobius transformations

Checking normality

Checking equicontinuity by Michael Becker

"Iteration is equicontinuous on the Fatou set and not on the Julia set". (Wolf Jung)

Michael Becker compares the distance of two close points under iteration on Riemann sphere.

This method can be used to draw not only Julia sets for polynomials ( where infinity is always superattracting fixed point) but it can be also applied to other functions ( maps), for which infinity is not an attracting fixed point.

using Marty's criterion by Wolf Jung

Wolf Jung is using "an alternative method of checking normality, which is based on Marty's criterion: |f'| / (1 + |f|^2) must be bounded for all iterates." It is implemented in mndlbrot::marty function ( see src code of program Mandel ver 5.3 ). It uses one point of dynamic plane.

Inverse Iteration Method (IIM/J) : Julia set

Inverse iteration of repellor for drawing Julia set

Complex potential - Boettcher coordinate

Exterior or complement of filled-in Julia set is :

• a basin of attraction of infinity ( superattracting fixed point)
• one of components of the Fatou set

$A_{f_c}(\infty) = K(f_c)^C = F_\infty$

It can be analysed using

• escape time (simple but gives only radial values = escape time ) LSM/J,
• distance estimation ( more advanced, continuus, but gives only radial values = distance ) DEM/J
• Boettcher coordinate or complex potential ( the best , gives  :
  • radial values ( real potential ) LogPhi = CPM/J
  • angular values ( external angle ) ArgPhi

So both values can be used to color with 2D gradient.

First compute Boettcher coordinate $w\,$ using this formula

$w = \Phi_c(z) = z* \prod_{n=1}^\infty \left( \frac{f_c^n(z)}{f_c^n(z)-c}\right )^{1/2^n}\,$

It looks "simple", but square root of complex nuber gives two values so one have to choose one value.

LogPhi - Douady-Hubbard potential - real potential - radial component of complex potential

CPM/J

Potential of filled Julia set

Diagram of potential computed with 2 methods : simple and full

Note that potential inside Kc is zero so :

Pseudocode version :

if (LastIteration==IterationMax)
 then potential=0    /* inside  Filled-in Julia set */
 else potential= GiveLogPhi(z0,c,ER,nMax); /* outside */

It also removes potential error for log(0).

Full version

Math (full) notation :

$LogPhi_c(z) = ln |z| + \sum_{n=1}^\infty \frac{1}{2^n} ln | 1 + \frac{c}{(f_c^{n-1}(z))^2}|$

Maxima (full) function :

GiveLogPhi(z0,c,ER,nMax):=
block(
 [z:z0,
 logphi:log(cabs(z)),
 fac:1/2,
 n:0],
 while n<nMax and abs(z)<ER do
 (z:z*z+c,
 logphi:logphi+fac*log(cabs(1+c/(z*z))),
 n:n+1
 ),
 return(float(logphi))
)$

Simplified version

The escape rate function of a polynomial f is defined by :

$G_f(z) = \lim_{n\rightarrow \infty}\frac{1}{2^n} log^+ |f^n(z)| \,$

where :

log ⁺ = max(log,0)

"The function Gp is continous on C and harmonic on the complement of the Julia set. It vanishes identically on K(f) and as it has a logarithmic pole at infinity, it is a it is the Green's function for C/ K(f)." ( Laura G. DeMarco)

Math simplified formula :

$SLogPhi_c(z) = \frac{log(f^n(z))}{ 2^n} \,$

Maxima function :

GiveSLogPhi(z0,c,e_r,i_max):=
 block(
  [z:z0,
   logphi,
   fac:1/2,
   i:0
  ],
  while i<i_max and cabs(z)<e_r do
  (z:z*z+c,
    fac:fac/2,
    i:i+1
  ),
 logphi:fac*log(cabs(z)),
 return(float(logphi))
)$

If you don't check if orbit is not bounded ( escapes, bailout test) then use this Maxima function :

GiveSLogPhi(z0,c,e_r,i_max):=
block(
 [z:z0, logphi, fac:1/2, i:0],
  while i<i_max and cabs(z)<e_r do
   (z:z*z+c,
    fac:fac/2,
    i:i+1 ),
   if i=i_max
    then logphi:0
    else logphi:fac*log(cabs(z)),
   float(logphi) 
)$

C version :

double jlogphi(double zx0, double zy0, double cx, double cy)
/* this function is based on function by W Jung http://mndynamics.com */
{ 
 int j; 
 double 
 zx=zx0,
 zy=zy0,
 s = 0.5, 
 zx2=zx*zx,
 zy2=zy*zy,
 t;
 for (j = 1; j < 400; j++)
 { s *= 0.5; 
  zy = 2 * zx * zy + cy;
  zx = zx2 - zy2 + cx; 
  zx2 = zx*zx; 
  zy2 = zy*zy;
  t = fabs(zx2 + zy2); // abs(z)
  if ( t > 1e24) break; 
 } 
return s*log2(t);  // log(zn)* 2^(-n)
}//jlogphi

Euler version by R. Grothmann ( with small change : from z^2-c to z^2+c) :

function iter (z,c,n=100) ...

h=z;
loop 1 to n;
h=h^2+c;
if totalmax(abs(h))>1e20; m=#; break; endif;
end;
return {h,m};
endfunction

x=-2:0.05:2; y=x'; z=x+I*y;
{w,n}=iter(z,c);
wr=max(0,log(abs(w)))/2^n;

Level Sets of potential = pLSM/J

Here is Delphi function which gives level of potential :

Function GiveLevelOfPotential(potential:extended):integer;
 var r:extended;
 begin
    r:= log2(abs(potential));
    result:=ceil(r);
 end;

ArgPhi - External angle - angular component of complex potential

One can start with binary decomposition of basin of attraction of infinity.

The second step can be using $\Psi_c\,$

polar coordinate system and Psi_c for c=-2

period detection

How to find period of external angle measured in turns under doubling map :

Here is Common Lisp code :

(defun give-period (ratio-angle)
  "gives period of angle in turns (ratio) under doubling map"
  (let* ((n (numerator ratio-angle))
         (d (denominator ratio-angle))
         (temp n)) ; temporary numerator
 
    (loop for p from 1 to 100 do 
          (setq temp  (mod (* temp 2) d)) ; (2 x n) modulo d = doubling)
          when ( or (= temp n) (= temp 0)) return p )))

Maxima CAS code :

doubling_map(n,d):=mod(2*n,d);

/* catch-throw version by Stavros Macrakis, works */
GivePeriodOfAngle(n0,d):=
catch(
      block([ni:n0],
          for i thru 200 do if (ni:doubling_map(ni,d))=n0 then throw(i),
          0 ) )$


/* go-loop version, works */
GiveP(n0,d):=block(
[ni:n0,i:0],
block(
  loop,
  ni:doubling_map(ni,d),
  i:i+1,
  if i<100 and not (n0=ni) then go(loop)
),
if (n0=ni) 
        then i 
        else 0
);

/* Barton Willis while version without for loop , works */
GivePeriod(n0,d):=block([ni : n0,k : 1],
  while (ni : doubling_map(ni,d)) # n0 and k < 100 do (
    k : k + 1),
  if k = 100 then 0 else k)$

Computing external angle

External angle (argument) is argument of Boettcher coordinate $w\,$

$arg_c(z) = arg(w) = arg(\Phi_c(z)) \,$

Because Boettcher coordinate is a product of complex numbers so argument of product is :

$arg_c(z) = arg(z) + \sum_{n=1}^\infty \left( \frac{1}{2^n}*arg \left( \frac{f_c^n(z)}{f_c^n(z)-c} \right ) \right ) \,$

Constructing the spine of filled Julia set

Algorithm for constructiong the spine is described by A. Douady

• join $- \beta\,$ and $\beta\,$,
• (to do )

Drawing dynamic external ray

Field lines in in the Fatou domain

Explanation by Gert Buschmann

backwards iteration

This method have been used by several people and proved by Thierry Bousch.

Code in c++ by Wolf Jung can be found in procedure QmnPlane::backray() in file qmnplane.cpp ( see source code of program mandel version 5.3 ).

• Ray for periodic angle ( simplest case )

It will be explained by an example :

First choose external angle $t \,$ (in turns). External angle for periodic ray is a rational number.

$t = \frac{1}{3}\,$

Compute period of external angle under doubling map.

Because "1/3 doubled gives 2/3 and 2/3 doubled gives 4/3, which is congruent to 1/3"

$1 \equiv 2 \pmod 3.\,$

or

$\frac{2}{3} \equiv \frac{1}{3} \pmod 1.\,$

so external angle $t = \frac{1}{3}\,$ has period 2 under doubling map.

Start with 2 points near infinity (in conjugate plane):

on ray 1/3 is a point $w = R*e^{2\pi i/3}\,$

on ray 2/3 is a point $w = R*e^{4 \pi i /3}\,$.

Near infinity $z = w \,$ so one can swith to dynamical plane ( Boettcher conjugation )

Backward iteration (with proper chose from two possibilities) of point on ray 1/3 goes to ray 2/3, back to 1/3 and so on.

Backward iteration (with proper chose from two possibilities) of point on ray 2/3 goes to ray 1/3, back to 2/3 and so on.

One has to divide set of points into 2 subsets ( 2 rays). Draw one of these 2 sets and join the points. It will be an approximation of ray.

• Ray for preperiodic angle ( to do )

Drawing dynamic external ray using inverse Boettcher map by Curtis McMullen

Julia set with external rays using McMullen method

This method is based on C program by Curtis McMullen and its Pascal version by Matjaz Erat

It consist of 3 big steps :

• compute some w-points of external ray of circle for angle $\theta \,$ and various radii (rasterisation)

$w_0 = R * e^{2\pi i t} \,$ where $1 < R < \infty \,$

• map w-points to z-point using inverse Boettcher map $\Psi_c = \Phi_{c}^{-1} \,$

$z_0 = \Psi_c (w_0) = f_c^{-n}(w_0^{2^n}) \,$

• draw z-points ( and connect them using lines )

First and last steps are easy, but second is not so needs more explanation.

Rasterisation

For given external ray in $f_0\,$ plane each point of ray has :

• constant value $t\,$ ( external angle in turns )
• variable radius $R\,$

so $w\,$ points of ray are parametrised by radius $R\,$ and can be computed using exponential form of complex numbers :

$w = F_t(R) = R * e^{2*\Pi*i*t}\,$

One can go along ray using linear scale :

t:1/3; /* example value */
R_Max:4;
R_Min:1.1;
for R:R_Max step -0.5 thru R_Min do w:R*exp(2*%pi*%i*t);
/* Maxima allows non-integer values in for statement */ 

It gives some w points with equal distance between them.

Another method is to use nonlinera scale.

To do it we introduce floating point exponent $r\,$ such that :

$R = 2^r\,$

and

$w = G_t(r) = 2^r * e^{2*\Pi*i*t}\,$

To compute some w points of external ray in $f_0\,$ plane for angle $t\,$ use such Maxima code :

t:1/3; /* external angle in turns */
/* range for computing  R ; as r tends to 0 R tends to 1 */
rMax:2; /* so Rmax=2^2=4 /
rMin:0.1; /* rMin > 0   */
caution:0.93; /* positive number < 1 ; r:r*caution  gives smaller r */
r:rMax;
unless r<rMin do
( 
 r:r*caution, /* new smaller r */ 
 R:2^r, /* new smaller R */
 w:R*exp(2*%pi*%i*t) /* new point w in f0 plane */
);

In this method distance between points is not equal but inversely proportional to distance to boundary of filled Julia set.

It is good because here ray has greater curvature so curve will be more smooth.

Mapping

Mapping points from $f_0\,$-plane to $f_c\,$-plane consist of 4 minor steps :

• forward iteration in $f_0\,$ plane

$w_n = w_0^{2^n} = R^{2^n}*e^{2\pi i 2^n t}\,$

until $w_n\,$ is near infinity

$abs(w_n) > R_{limit} \,$

• switching plane ( from $f_0\,$ to $f_c\,$)

$z_n =w_n \,$

( because here, near infinity : $w = z \,$)

• backward iteration in $f_c\,$ plane the same ( $n \,$) number of iterations
• last point $z_0 \,$ is on our external ray

1,2 and 4 minor steps are easy. Third is not.

$z_{n-1} = f_c^{-1} (z_n) = \sqrt{z_n - c}\,$

Backward iteration uses square root of complex number. It is 2-valued functions so backward iteration gives binary tree.

One can't choose good path in such tree without extre informations. To solve it we will use 2 things :

• equicontinuity of basin of attraction of infinity
• conjugacy between $f_0\,$ and $f_c\,$ planes

Equicontinuity of basin of attraction of infinity

Basin of attraction of infinity ( complement of filled-in Julia set) contains all points which tends to infinity under forward iteration.

$A_{f_c}(\infty) \ \overset{\underset{\mathrm{def}}{}}{=} \ \{ z \in \mathbb{C}  : f^{(k)} _c (z) \to \infty\ as\ k \to \infty \}.$

Infinity is superattracting fixed point and orbits of all points have similar behaviour. In other words orbits of 2 points are assumed to stay close if they are close at the beginning.

It is equicontinuity ( compare with normality).

In $f_c\,$ plane one can use forward orbit of previous point of ray for computing backward orbit of next point.

Detailed version of algorithm

• compute first point of ray (start near infinity ang go toward Julia set )

$w = 2^r * e^{2\Pi i t}\,$ where $r = r_{max}\,$

here one can easily switch planes :

$z = w\,$

It is our first z-point of ray.

• compute next z-point of ray
  • compute next w-point of ray for $r = r * caution \,$
  • compute forward iteration of 2 points : previous z-point and actual w-point. Save z-orbit and last w-point
  • switch planes and use last w-point as a starting point :$z_n = w_{last}\,$
  • backward iteration of new $z_n\,$ toward new $z_0\,$ using forward orbit of previous z point
  • $z_0\,$ is our next z point of our ray
• and so on ( next points ) until $r < r_{Min}\,$

Maxima CAS src code

/* gives a list of z-points of external ray for angle t in turns and coefficient c */
GiveRay(t,c):=
block(
 [r],
 /* range for drawing  R=2^r ; as r tends to 0 R tends to 1 */
 rMin:1E-20, /* 1E-4;  rMin > 0  ; if rMin=0 then program has infinity loop !!!!! */
 rMax:2,  
 caution:0.9330329915368074, /* r:r*caution ; it gives smaller r */
 /* upper limit for iteration */
 R_max:300,
 /* */
 zz:[], /* array for z points of ray in fc plane */
 /*  some w-points of external ray in f0 plane  */
 r:rMax,
 while 2^r<R_max do r:2*r, /* find point w on ray near infinity (R>=R_max) in f0 plane */
 R:2^r,
 w:rectform(ev(R*exp(2*%pi*%i*t))),
 z:w, /* near infinity z=w */
 zz:cons(z,zz), 
 unless r<rMin do
 ( /* new smaller R */
  r:r*caution,  
  R:2^r,
  /* */
  w:rectform(ev(R*exp(2*%pi*%i*t))),
  /* */
  last_z:z,
  z:Psi_n(r,t,last_z,R_max), /* z=Psi_n(w) */
  zz:cons(z,zz)
 ),
 return(zz)
)$

Lamination of Julia set

Lamination of rabbit Julia set

• Douady’s pinched disks
• Thurston’s invariant laminations
• Milnor’s orbit portraits

Here is long description

BSM/J

This algorithm is used when dynamical plane consist of two of more basins of attraction. For example for c=0.

It is not appropiate when interior of filled Julia set is empty, for example for c=i.

Description of algorithm :

• for every pixel of dynamical plane z do :
  • compute 4 corners ( vertices) of pixel z_lt,z_rt,z_rb,z_lb ( where lt denotes left top, rb denotes right bottom, ... )
  • check to which basin corner belongs ( standard escape time and bailout test )
  • if corners do not belong to the same basin mark it as Julia set

Examples of code

• program in Pascal
• via convolution with a kernel

DEM/J

External distance estimation

• 

  c=-0.74543+0.11301*i

• 

  c= -0.75+0.11

• 

  c = 0.255

• 

  c=-0.1+0.651

" For distance estimate it has been proved that the computed value differs from the true distance at most by a factor of 4. " (Wolf Jung)

Math formula :

$distance(z_0, K_c, n) = 2 * |z_n| * \frac{log|z_n| }{ |z'_n| }\,$

where :

$z'_n\,$ is first derivative with respect to c.

This derivative can be found by iteration starting with

$z'_0 = 1 \,$

and then :

$z'_n= 2*z_{n-1}*z'_{n-1}\,$

How to use distance

One can use distance for colouring  :

• only Julia set ( boundary of filled Julia set)
• boundary and exterior of filled Julia set.

Here is first example :

if (LastIteration==IterationMax) 
    then { /*  interior of Julia set, distance = 0 , color black */ }
    else /* exterior or boundary of Filled-in Julia set  */
         {  double distance=give_distance(Z0,C,IterationMax);
            if (distance<distanceMax)
                then { /*  Julia set : color = white */ }
                else  { /*  exterior of Julia set : color = black */}
         }

Here is second example

if (LastIteration==IterationMax) or distance < distanceMax then ... // interior by ETM/J and boundary by DEM/J
else .... // exterior by real escape time

Zoom

DistanceMax is smaller than pixel size. During zooming pixel size is decreasing and DistanceMax should also be decreased to obtain good picture. It can be made by using formula :

$DistanceMax = \frac{PixelSize}{n} \,$

where $n \le 1 \,$

One can start with n=1 and increase n if picture is not good. Check also iMax !!

DistanceMax may also be proportional to zoom factor $mag\,$ :

$DistanceMax = \sqrt{\frac {thick}{1000 * mag}}$

where thick is image width ( in world units) and mag is a zoom factor.

Examples of code

In cpp see mndlbrot::dist from mndlbrot.cpp in src code of program mandel ver 5.3 by Wolf Jung.

C function :

/*based on function  mndlbrot::dist  from  mndlbrot.cpp
 from program mandel by Wolf Jung (GNU GPL )
 http://www.mndynamics.com/indexp.html  */
double jdist(double Zx, double Zy, double Cx, double Cy ,  int iter_max)
 { 
 int i;
 double x = Zx, /* Z = x+y*i */
         y = Zy, 
         /* Zp = xp+yp*1 = 1  */
         xp = 1, 
         yp = 0, 
         /* temporary */
         nz,  
         nzp,
         /* a = abs(z) */
         a; 
 for (i = 1; i <= iter_max; i++)
  { /* first derivative   zp = 2*z*zp  = xp + yp*i; */
    nz = 2*(x*xp - y*yp) ; 
    yp = 2*(x*yp + y*xp); 
    xp = nz;
    /* z = z*z + c = x+y*i */
    nz = x*x - y*y + Cx; 
    y = 2*x*y + Cy; 
    x = nz; 
    /* */
    nz = x*x + y*y; 
    nzp = xp*xp + yp*yp;
    if (nzp > 1e60 || nz > 1e60) break;
  }
 a=sqrt(nz);
 /* distance = 2 * |Zn| * log|Zn| / |dZn| */
 return 2* a*log(a)/sqrt(nzp); 
 }

Delphi function :

function Give_eDistance(zx0,zy0,cx,cy,ER2:extended;iMax:integer):extended;
 
var zx,zy ,  // z=zx+zy*i
    dx,dy,  //d=dx+dy*i  derivative :  d(n+1)=  2 * zn * dn
    zx_temp,
    dx_temp,
    z2,  //
    d2, //
    a // abs(d2)
     :extended;
    i:integer;
begin
  //initial values
  // d0=1
  dx:=1;
  dy:=0;
  //
  zx:=zx0;
  zy:=zy0;
  // to remove warning : variables may be not initialised ?
  z2:=0;
  d2:=0;
 
  for i := 0 to iMax - 1 do
    begin
      // first derivative   d(n+1) = 2*zn*dn  = dx + dy*i;
      dx_temp := 2*(zx*dx - zy*dy) ;
      dy := 2*(zx*dy + zy*dx);
      dx := dx_temp;
      // z = z*z + c = zx+zy*i
      zx_temp := zx*zx - zy*zy + Cx;
      zy := 2*zx*zy + Cy;
      zx := zx_temp;
      //
      z2:=zx*zx+zy*zy;
      d2:=dx*dx+dy*dy;
      if ((z2>1e60) or (d2 > 1e60)) then  break;
 
    end; // for i
   if (d2 < 0.01) or (z2 < 0.1)  // when do not use escape time
    then  result := 10.0
    else
      begin
        a:=sqrt(z2);
        // distance = 2 * |Zn| * log|Zn| / |dZn|
        result := 2* a*log10(a)/sqrt(d2);
      end;
 
end;  //  function Give_eDistance

Matlab code by Jonas Lundgren

function D = jdist(x0,y0,c,iter,D)
%JDIST Estimate distances to Julia set by potential function
% by Jonas Lundgren http://www.mathworks.ch/matlabcentral/fileexchange/27749-julia-sets
% Code covered by the BSD License http://www.mathworks.ch/matlabcentral/fileexchange/view_license?file_info_id=27749
 
 
% Escape radius^2
R2 = 100^2;
 
% Parameters
N = numel(x0);
M = numel(y0);
cx = real(c);
cy = imag(c);
iter = round(1000*iter);
 
% Create waitbar
h = waitbar(0,'Please wait...','name','Julia Distance Estimation');
t1 = 1;
 
% Loop over pixels
for k = 1:N/2
    x0k = x0(k);
    for j = 1:M
        % Update distance?
        if D(j,k) == 0
            % Start values
            n = 0;
            x = x0k;
            y = y0(j);
            b2 = 1;                 % |dz0/dz0|^2
            a2 = x*x + y*y;         % |z0|^2
            % Iterate zn = zm^2 + c, m = n-1
            while n < iter && a2 <= R2
                n = n + 1;
                yn = 2*x*y + cy;
                x = x*x - y*y + cx;
                y = yn;
                b2 = 4*a2*b2;       % |dzn/dz0|^2
                a2 = x*x + y*y;     % |zn|^2
            end
            % Distance estimate
            if n < iter
                % log(|zn|)*|zn|/|dzn/dz0|
                D(j,k) = 0.5*log(a2)*sqrt(a2/b2);
            end
        end
    end
    % Lap time
    t = toc;
    % Update waitbar
    if t >= t1
        str = sprintf('%0.0f%% done in %0.0f sec',200*k/N,t);
        waitbar(2*k/N,h,str)
        t1 = t1 + 1;
    end
end
 
% Close waitbar
close(h)

Maxima function :

GiveExtDistance(z0,c,e_r,i_max):=
/* needs z in exterior of Kc */
block(
[z:z0,
dz:1,
cabsz:cabs(z),
cabsdz:1, /* overflow limit */
i:0],
while  cabsdz < 10000000 and  i<i_max
 do 
  (
   dz:2*z*dz,
   z:z*z + c,
   cabsdz:cabs(dz),
   i:i+1
  ),
 cabsz:cabs(z), 
 return(2*cabsz*log(cabsz)/cabsdz)
)$

Cauchy Convergence Algorithm (CCA)

This algorithm is described by User:Georg-Johann. Here is also Matemathics code by Paul Nylander

Koenigs coordinate

Koenigs coordinate are used in the basin of attraction of finite attracting (not superattracting) point (cycle),

Optimisation

Symmetry

Julia sets can have many symmetries

Quadratic Julia set has allways rotational symmetry ( 180 degrees) :

colour(x,y) = colour(-x,-y)

when c is on real axis ( cy = 0) Julia set is also reflection symmetric :

colour(x,y) = colour(x,-y)

Algorithm :

• compute half image
• rotate and add the other half
• write image to file

Color

• Visualising Julia sets by Georg-Johann
• Combined Methods of Depicting Julia Sets and Parameter Planes by Chris King
• On Fractal Coloring Techniques by Jussi Harkonen Master's Thesis, Department of Mathematics, Åbo Akademi University, Turku, 2007, 61 pages. The thesis was carried out under the supervision of Professor Goran Hognas
• Technical Info - Colorizing by Michael Condron
• Colors/Color_gradient
• Technicolor Julias by Shawn Hargreaves

Sets

Target set

For infinite attractor

Target set $T\,$ is an arbitrary set on dynamical plane containing infinity and not containing points of Filled-in Fatou sets.

For escape time algorithms target set determines the shape of level sets and curves. It does not do it for other methods.

Exterior of circle

This is typical target set. It is exterior of circle with center at origin $z = 0 \,$ and radius =ER :

$T_{ER}=\{z:abs(z) > ER \} \,$

Radius is named escape radius ( ER ) or bailout value. Radius should be greater than 2.

Exterior of square

Here target set is exterior of square of side length $s\,$ centered at origin

$T_s=\{z: abs(re(z)) > s ~~\mbox{or}~~ abs(im(z))>s \} \,$

Julia sets

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

• when Julia set is a set of points that do not escape to infinity under iteration of the quadratic map ( = filled Julia set has no interior = dendrt)
  • IIM/J
  • DEM/J
  • checking normality
• when Julia set is a boundary between 2 basin of attraction ( = filled Julia set has no empty interior) :
  • boundary scaning
  • edge detection

Fatou set

Interior of filled Julia set can be coloured :

• speed of attraction ( integer value = the number of iterations used to guess if a point is in the set ) which is coverted to colour ( or shade of gray )
  • attracting time ( sth like escape time but checks if (abs(z-attractor)<Attracting_radius
• discrete veolocity in Siegel disc case

Periodic points

Finding equations for periodic points

Standard equations

Equation for periodic point for period $p \,$ is :

$\ F_p(z,c) = 0 \,$

so in Maxima :

for i:1 thru 4 step 1 do disp(i,expand(Fn(i,z,c)=0));
1
z^2-z+c=0
2
z^4+2*c*z^2-z+c^2+c=0
3
z^8+4*c*z^6+6*c^2*z^4+2*c*z^4+4*c^3*z^2+4*c^2*z^2-z+c^4+2*c^3+c^2+c=0
4
z^16+8*c*z^14+28*c^2*z^12+4*c*z^12+56*c^3*z^10+24*c^2*z^10+70*c^4*z^8+60*c^3*z^8+6*c^2*z^8+2*c*z^8+56*c^5*z^6+80*c^4*z^6+
24*c^3*z^6+8*c^2*z^6+28*c^6*z^4+60*c^5*z^4+36*c^4*z^4+16*c^3*z^4+4*c^2*z^4+8*c^7*z^2+24*c^6*z^2+24*c^5*z^2+16*c^4*z^2+
8*c^3*z^2- z+c^8+4*c^7+6*c^6+6*c^5+5*c^4+2*c^3+c^2+c=0

Degree of standard equations is $2^p \,$

These equations give periodic points for period p and its divisors. It is because these polynomials are product of lower period polynomials $G_p\,$:

for i:1 thru 4 step 1 do disp(i,factor(Fz(i,z,c)));
1
z^2-z+c
2
(z^2-z+c)*(z^2+z+c+1)
3
(z^2-z+c)*(z^6+z^5+3*c*z^4+z^4+2*c*z^3+z^3+3*c^2*z^2+3*c*z^2+z^2+c^2*z+2*c*z+z+c^3+2*c^2+c+1)
4
(z^2-z+c)*(z^2+z+c+1)(z^12+6*c*z^10+z^9+15*c^2*z^8+3*c*z^8+4*c*z^7+20*c^3*z^6+12*c^2*z^6+z^6+6*c^2*z^5+2*c*z^5+15*c^4*z^4+18*c^3*z^4+3*
c^2*z^4+4*c*z^4+4*c^3*z^3+4*c^2*z^3+z^3+6*c^5*z^2+12*c^4*z^2+6*c^3*z^2+5*c^2*z^2+c*z^2+c^4*z+2*c^3*z+c^2*z+2*c*z+c^6+
3*c^5+3*c^4+3*c^3+2*c^2+1)

so

$\ F_1(z,c) = z^2-z+c = G_1 \,$

$\ F_2(z,c) = (z^2-z+c)*(z^2+z+c+1) = G_1 * G_2\,$

$\ F_3(z,c) = (z^2-z+c)*(z^6+z^5+3*c*z^4+z^4+2*c*z^3+z^3+3*c^2*z^2+3*c*z^2+z^2+c^2*z+2*c*z+z+c^3+2*c^2+c+1) = G_1 * G_3\,$

$\ F_4(z,c) = G_1 * G_2 * G_4 \,$

$\ F_4(z,c) = G_1 * G_5 \,$

Standard equations can be reduced to equations $G_p\,$ giving periodic points for period p without its divisors.

Reduced equations

Reduced equation is

$G_p(z,c) = 0 \,$

so in Maxima gives :

for i:1 thru 5 step 1 do disp(i,expand(G[i,z,c]=0));
1
z^2-z+c=0
2
z^2+z+c+1=0
3
z^6+z^5+3*c*z^4+z^4+2*c*z^3+z^3+3*c^2*z^2+3*c*z^2+z^2+c^2*z+2*c*z+z+c^3+2*c^2+c+1=0
4
z^12+6*c*z^10+z^9+15*c^2*z^8+3*c*z^8+4*c*z^7+20*c^3*z^6+12*c^2*z^6+z^6+6*c^2*z^5+2*c*z^5+15*c^4*z^4+18*c^3*z^4+
3*c^2*z^4+4*c*z^4+4*c^3*z^3+4*c^2*z^3+z^3+6*c^5*z^2+12*c^4*z^2+
6*c^3*z^2+5*c^2*z^2+c*z^2+c^4*z+2*c^3*z+c^2*z+2*c*z+c^6+3*c^5+3*c^4+3*c^3+2*c^2+1=0

Computing periodic points in Maxima

Define c value

(%i1) c:%i;
(%o1) %i

Fixed points ( period = 1 )

Compute fixed points

(%i2) p:float(rectform(solve([z*z+c=z],[z])));
(%o2) [z=0.62481053384383*%i-0.30024259022012,z=1.30024259022012-0.62481053384383*%i]

Find which is repelling :

abs(2*rhs(p[1]));

Show that sum of fixed points is 1

(%i10) p:solve([z*z+c=z], [z]);
(%o10) [z=-(sqrt(1-4*%i)-1)/2,z=(sqrt(1-4*%i)+1)/2]
(%i14) s:radcan(rhs(p[1]+p[2]));
(%o14) 1

Draw points :

(%i15) xx:makelist (realpart(rhs(p[i])), i, 1, length(p));
(%o15) [-0.30024259022012,1.30024259022012]
(%i16) yy:makelist (imagpart(rhs(p[i])), i, 1, length(p));
(%o16) [0.62481053384383,-0.62481053384383]
(%i18) plot2d([discrete,xx,yy], [style, points]);

One can add point z=1/2 to the lists:

(%i31) xx:cons(1/2,xx);
(%o31) [1/2,-0.30024259022012,1.30024259022012]
(%i34) yy:cons(0,yy);
(%o34) [0,0.62481053384383,-0.62481053384383]
(%i35) plot2d([discrete,xx,yy], [style, points]);

period 2 points

(%i2) solve([(z*z+c)^2+c=z], [z]);
(%o2) [z=-(sqrt(-4*c-3)+1)/2,z=(sqrt(-4*c-3)-1)/2,z=-(sqrt(1-4*c)-1)/2,z=(sqrt(1-4*c)+1)/2]

Show that z1+z2 = -1

(%i4) s:radcan(rhs(p[1]+p[2]));
(%o4) -1

Show that z2+z3=1

(%i6) s:radcan(rhs(p[3]+p[4]));
(%o6) 1

Finding period of orbit

critical orbit with 3-period cycle

Drawing Julia set and critical orbit. Computing period

"methods based on interval arithmetic when implemented properly are capable of finding all period-n cycles for considerable large n." (ZBIGNIEW GALIAS )

Models of Julia set

• lamination of Julia set
• graph of Julia set
  • Schreier graphs
• spine of Julia set
• Hubbard trees serve as a model for Julia sets of post-critically finite polynomials

Video

One can make videos using :

• zoom into dynamic plane
• changing parametr c along path inside parameter plane
• changing coloring scheme ( for example color cycling )

More tutorials and code

• in Java see Evgeny Demidov
• in C see :
  • art gallery of Linas Vepstas
  • Fractint
  • Gfract: multi-threaded fractal exploration program by Osku Salerma
  • mdz : Mandelbrot Deep Zoom open source software for Linux by James W. Morris, based on Gfract
• in C++ see Wolf Jung page,
• in Gnuplot see Tutorial by T.Kawano
• in Lisp for Maxima see Dynamics by Jaime E. Villate
• in Mathemathica see :
  • Inverse Iteration Algorithms for Julia Sets by Mark McClure
  • Robert M. Dickau code

References

• Drakopoulos V., Comparing rendering methods for Julia sets, Journal of WSCG 10 (2002), 155–161

• tree with dynamics by Nathaniel D. Emerson