From Wikibooks, the open-content textbooks collection

Current revision (unreviewed)

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

[edit] Algorithms

[edit] Koenigs coordinate

Koenigs^[1] coordinate are used in the basin of attraction of finite attracting (not superattracting) point (cycle),^[2]^[3]

[edit] Escape time

[edit] Basin of attraction to infinity = exterior of filled-in Julia set

Here one computes forward iterations of a complex point Z_n:

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

Here is function which computes the last iteration, that is the first iteration that 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

^[4]

Lisp version

((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 :

GiveLastIteration(z,c):=
 block([i:0],
  while abs(z)<ER and i<iMax
     do (z:z*z + c,i:i+1),
  i)$

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;

[edit] Boolean Escape time

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

Algorithm: for every point z of dynamical plane (z-plane) compute iteration number ( last iteration) for which magnitude of z is greater then 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 .

[edit] 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;

[edit] Normalized iteration count (real escape time)

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-log(log(cabs(z))/log(2)),
 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);

[edit] Level Curves of escape time Method = eLCM/J

Level Curves of Escape time. Image and C code

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

edge detection of Level sets. Algorithm is based on paper by M. Romera et al.^[5]
drawing lemniscates = curves $L_n=\{z: abs(z_n)=ER \}\,$ , see explanation and source code
method described by Harold V. McIntosh^[6]

/* 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]);

[edit] 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

Interior of filled-in Julia set for C = 0.5*i. Image and source code

[edit] 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^[7]

Components of Fatou set. Image and source code

[edit] Decomposition of level sets

[edit] Binary decomposition

Fatou set for F(z)=z*z made with binary decomposition Image and source code

Binary decomposition of basin of attraction of infinity (algorithm in pseudocode) :

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=Re(Z) */
           then color=BLACK; 
           else color=WHITE;

[edit] Checking normality

[edit] Checking equicontinuity by Michael Becker

"Iteration is equicontinuous on the Fatou set and not on the Julia set". (Wolf Jung) ^[8]^[9]

Michael Becker compares the distance of two close points under iteration on Riemann sphere.^[10]

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

[edit] 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.

[edit] Inverse Iteration Method (IIM/J) : Julia set

Julia set for C = i. Image and source code

Julia set made with MIIM. Image and maxima source code

distribution of points in simple IIM/J

In escape time one computes forward iteration of point z.

In IIM/J one computes:

repelling fixed point of complex quadratic polynomial $Z_0=\beta_c \,$
preimages of $Z_0\,$ by inverse iterations

$Z_{n-1} = \sqrt{Z_n - C}$

Because square root is multivalued function then each $Z_{n}\,$ has two preimages $Z_{n-1}\,$ . Thus inverse iteration creates binary tree.

Here is Maxima code for simple IIM:

finverseplus(z,c):=sqrt(z-c);
finverseminus(z,c):=-sqrt(z-c);
c:%i;    /*       define c value  */
iMax:5000;  /* maximal number of reversed iterations */
fixed:float(rectform(solve([z*z+c=z],[z])));  /* compute fixed points of f(z,c) :   z=f(z,c)   */
if (abs(2*rhs(fixed[1]))<1)  /* Find which is repelling  */
then block(beta:rhs(fixed[1]),alfa:rhs(fixed[2])) 
else block(alfa:rhs(fixed[1]),beta:rhs(fixed[2]));
 z:beta; /* choose repeller as a starting point */
/* make 2 list of points and copy beta to  to lists */ 
xx:makelist (realpart(beta), i, 1, 1); /* list of re(z) */
yy:makelist (imagpart(beta), i, 1, 1); /* list of im(z) */ 
for i:1 thru iMax  step 1 do  /* reversed iteration of beta */
block
(if random(1.0)>0.5
  then z:finverseplus(z,c)
  else z:finverseminus(z,c),
 xx:cons(realpart(z),xx), /* save values to draw it later */
 yy:cons(imagpart(z),yy)
);
plot2d([discrete,xx,yy],[style,[points,1,0,1]]); /* draw reversed orbit of beta */

Compare it with C code.

[edit] Variants of IIM : MIIM/J

random choose one of two roots ( up to chosen level )= random walk through the tree.
- both roots with the same probability
- more often one then other root
draw some rare paths in binary tree = MIIM
draw all roots ( up to chosen level ). There are some posible way to go through all the tree. It is implemented in Fractint^[12].

[edit] 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 Wolf Jung formula ^[13]

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

[edit] LogPhi - Douady-Hubbard potential - real potential - radial component of complex potential

Potential of filled Julia set

Diagram of potential computed with 2 methods : simple and full

Note that potential inside Kc is zero so (in pseudocode ) :

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

Math (full) version : ^[14]

$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) version :

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

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 dont 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 */
{ 
 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

[edit] CPM/J

[edit] 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;

[edit] Level Curves of potential = equipotential lines = pLCM/J

[edit] 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

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 ) \,$

[edit] Constructing the spine of filled Julia set

Algorithm for constructiong the spine is described by A. Douady^[15]

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

[edit] Drawing dynamic external ray

[edit] backwards iteration

This method have been used by several people and proved by Thierry Bousch. ^[16]

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 ). ^[17]

Periodic ray ( 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" ^[18]

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

or

$\frac{2}{3} \equiv \frac{1}{3} \pmod 1.\,$ ^[19]

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 $w = R*e^{2\pi i/3}\,$ and on ray 2/3 $w = R*e^{4 \pi i /3}\,$ , here $z = w \,$ so one can swith to dynamical plane ( Boettcher conjugation )

Backward iteration (with proper chose from two possibilities)^[20] 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 points of backward orbit between 2 rays and can draw this 2 rays.

Preperiodic ray

[edit] 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^[21] and its Pascal version by Matjaz Erat^[22]

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 conect them using lines )

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

[edit] 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.

[edit] 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

[edit] =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.

[edit] 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}\,$

[edit] 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)
)$

[edit] BSM/J

[edit] DEM/J

[edit] Internal distance estimation

[edit] External distance estimation

Julia set : image with C source code

" 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}\,$

Algorithm for computing distance from point z ( in exterior ) to nearest point on the boundary of Julia set. Here is pseudocode:

if (LastIteration==IterationMax) 
    then { /*  interior of Julia set, distance = 0 , color black */ }
    else /* exterior 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 */}
         }

DistanceMax is smaller then pixel size. One can start with DistanceMax=pixel_size and if the result is not good change DistanceMax=pixel_size/n where n is number greater then 1.

Examples of code :

in cpp see mndlbrot::dist from mndlbrot.cpp in src code of program mandel ver 5.3 by Wolf Jung.
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)
)$

[edit] 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

[edit] Sets

[edit] 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 )^[23]

IIM/J
DEM/J

[edit] Fatou set

[edit] Periodic points

[edit] Finding equations for periodic points

[edit] 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 it's divisors.

[edit] 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

[edit] Computing periodic points in Maxima

Define c value

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

[edit] 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]);

[edit] 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

[edit] Finding period of orbit

critical orbit with 3-period cycle

Drawing Julia set and critical orbit. Computing period

[edit] More tutorials and code

in Java see Evgeny Demidov
in C see :
- art gallery of Linas Vepstas
- Fractint
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

[edit] References

↑ Gabriel Koenigs biographie at The MacTutor History of Mathematics archive
↑ G. Koenigs, Recherches sur les intégrales de certaines équations fonctionnelles, Annales École Normale Supérieure, 1(3) (1884), Supplément, 3-41.
↑ Inigo Quilez images and tutuorial
↑ Program Mandel by Wolf Jung on GNU General Public License
↑ Drawing the Mandelbrot set by the method of escape lines. M. Romera et al.
↑ Julia Curves, Mandelbrot Set, Harold V. McIntosh.
↑ The fixed points and periodic orbits by E Demidov
↑ Alan F. Beardon, S. Axler, F.W. Gehring, K.A. Ribet : Iteration of Rational Functions: Complex Analytic Dynamical Systems. Springer, 2000; ISBN 0387951512, 9780387951515; page 49
↑ Joseph H. Silverman : The arithmetic of dynamical systems. Springer, 2007. ISBN 0387699031, 9780387699035; page 22
↑ Visualising Julia sets by Georg-Johann
↑ Julia sets by Michael Becker. See the metric d(z,w)
↑ Fractint documentation - Inverse Julias
↑ How to draw external rays by Wolf Jung
↑ The Beauty of Fractals, page 65
↑ A. Douady, “Algorithms for computing angles in the Mandelbrot set,” in Chaotic Dynamics and Fractals, M. Barnsley and S. G. Demko, Eds., vol. 2 of Notes and Reports in Mathematics in Science and Engineering, pp. 155–168, Academic Press, Atlanta, Ga, USA, 1986.
↑ Thierry Bousch : De combien tournent les rayons externes? Manuscrit non publié, 1995
↑ Program Mandel by Wolf Jung
↑ Explanation by Wolf Jung
↑ Modular arithmetic in wikipedia
↑ Square root of complex number gives 2 values so one has to choose only one. For details see Wolf Jung page
↑ c program by Curtis McMullen (quad.c in Julia.tar.gz)
↑ Quadratische Polynome by Matjaz Erat
↑ Mark Braverman : On efficient computation of parabolic Julia sets

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

[0] Gabriel Koenigs biographie at The MacTutor History of Mathematics archive

[1] G. Koenigs, Recherches sur les intégrales de certaines équations fonctionnelles, Annales École Normale Supérieure, 1(3) (1884), Supplément, 3-41.

[2] Inigo Quilez images and tutuorial

[3] Program Mandel by Wolf Jung on GNU General Public License

[4] Drawing the Mandelbrot set by the method of escape lines. M. Romera et al.

[5] Julia Curves, Mandelbrot Set, Harold V. McIntosh.

[6] The fixed points and periodic orbits by E Demidov

[7] Alan F. Beardon, S. Axler, F.W. Gehring, K.A. Ribet : Iteration of Rational Functions: Complex Analytic Dynamical Systems. Springer, 2000; ISBN 0387951512, 9780387951515; page 49

[8] Joseph H. Silverman : The arithmetic of dynamical systems. Springer, 2007. ISBN 0387699031, 9780387699035; page 22

[9] Visualising Julia sets by Georg-Johann

[10] Julia sets by Michael Becker. See the metric d(z,w)

[11] Fractint documentation - Inverse Julias

[12] How to draw external rays by Wolf Jung

[13] The Beauty of Fractals, page 65

[14] A. Douady, “Algorithms for computing angles in the Mandelbrot set,” in Chaotic Dynamics and Fractals, M. Barnsley and S. G. Demko, Eds., vol. 2 of Notes and Reports in Mathematics in Science and Engineering, pp. 155–168, Academic Press, Atlanta, Ga, USA, 1986.

[15] Thierry Bousch : De combien tournent les rayons externes? Manuscrit non publié, 1995

[16] Program Mandel by Wolf Jung

[17] Explanation by Wolf Jung

[18] Modular arithmetic in wikipedia

[19] Square root of complex number gives 2 values so one has to choose only one. For details see Wolf Jung page

[20] rogram by Curtis McMullen (quad.c in Julia.tar.gz)

[21] Quadratische Polynome by Matjaz Erat

[22] Mark Braverman : On efficient computation of parabolic Julia sets

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

Fractals/Iterations in the complex plane/Julia set