# 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[1]
• 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.

How to find:

• lowest optimal bailout values ( IterationMax) ? [2]

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

Here one computes forward iterations of a complex point Z0:

${\displaystyle 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). It can also be improved [4]

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 with complex type from GSL :[5]

#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <stdio.h>
// gcc -L/usr/lib -lgsl -lgslcblas -lm t.c
// function fc(z) = z*z+c

gsl_complex f(gsl_complex z, gsl_complex c) {
}

int main () {
gsl_complex c = gsl_complex_rect(0.123, 0.125);
gsl_complex z = gsl_complex_rect(0.0, 0.0);
int i;
for (i = 0; i < 10; i++) {
z = f(z, c);
double zx = GSL_REAL(z);
double zy = GSL_IMAG(z);
printf("Real: %f4 Imag: %f4\n", zx, zy);
}
return 0;
}


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
{
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) [7]

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

((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 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: 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;  Here is the c function which: • uses complex double type numbers • computes 8 bit color ( shades of gray) • checks both escape and attraction test unsigned char ComputeColorOfLSM(complex double z){ int nMax = 255; double cabsz; unsigned char iColor; int n; for (n=0; n < nMax; n++){ //forward iteration cabsz = cabs(z); if (cabsz > ER) break; // esacping if (cabsz< PixelWidth) break; // fails into finite attractor = interior z = z*z +c ; /* forward iteration : complex quadratic polynomial */ } iColor = 255 - 255.0 * ((double) n)/20; // nMax or lower walues in denominator return iColor; }   "if a 2-variable function z = f(x,y) has non-extremal critical points, i.e. it has saddle points, then it's best if the contour z heights are chosen so that the saddle points are on a contour, so that the crossing contours appear visually."Alan Ableson  How to choose parameters for which Level curves cross critical point ( and it's preimages ) ? • choose the parameter c such that it is on an escape line, then the critical value will be on an escape line as well • choose escape radius equal to n=th iteration of critical value // find such ER for LSM/J that level curves croses critical point and it's preimages ( only for disconnected Julia sets) double GiveER(int i_Max){ complex double z= 0.0; // criical point int i; ; // critical point escapes very fast here. Higher valus gives infinity for (i=0; i< i_Max; ++i ){ z=z*z +c; } return cabs(z); }  #### Normalized iteration count (real escape time or fractional iteration or Smooth Iteration Count Algorithm (SICA)) Math formula : ${\displaystyle \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);


description:

#### Level Curves of escape time Method = eLCM/J

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.[8]
• Sobel filter
• drawing lemniscates = curves ${\displaystyle L_{n}=\{z:abs(z_{n})=ER\}\,}$, see explanation and source code
• drawing circle ${\displaystyle L_{0}=\{z:abs(z)=ER\}\,}$ and its preimages. See this image, explanation and source code
• method described by Harold V. McIntosh[9]
/* 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]);


Density of level curves[10]

  "The spacing between level curves is a good way to estimate gradients: level curves that are close together represent areas of steeper descent/ascent." [11]


 "The density of the contour lines tells how steep is the slope of the terrain/function variation. When very close together it means f is varying rapidly (the elevation increase or decrease rapidly). When the curves are far from each other the variation is slower" [12]


### 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[13]
• use edge detection for drawing Julia set

#### Internal Level Sets

See :

• algorithm 0 of program Mandel by Wolf Jung

How to choose size of attracting trap

• petal in the parabolic case :
• for period 1 and 2 : radius of a circle with parabolic point on it's boundary
• for higher periods sectore of the circle with center at parabolic point
• radius of the circle with attractor as a center

such that level curves cross at critical point ?

/*
https://commons.wikimedia.org/wiki/File:Fatou_sets_for_Blaschke_fraction_f(z)_%3D_rho_*_z%5E2_*_(z-3)_over_(1-3z)_with_LCM_and_critical_orbit.png
find AR  and center of parabolic trap ( sepal) by simply iteration of critical point
uses global var AR and
*/
int  GiveTunedAR_andTrapCenter(const complex double zcr, const complex double zf, const int iter_Max, const int period){

complex double z = zcr; // initial point z0 = criical point

// iterate critical point
for (int iter=0; iter< iter_Max; ++iter ){

for (int p =0; p< period; ++p )
{z = f(z);} // forward iteration

}
// check distance between zn = f^n(zcr) and periodic point zp0

//fprintf(stdout, "critical point zn =  %.16f %+.16f*I \n", creal (z), cimag (z));
AR = cabs(z - zf)/2.0; // !!!!
trap_center =  (creal(z) + creal(zf))/2.0 + I*(cimag(zf) + cimag(z))/2.0; // midpoint between zf and z

return 0;

}


For weakly attracting :

// compute radius of circle around finite attractor which is independent of the image size ( iWidth/2000.0 )
// input k is a number of pixels ( in case of iWidth = 2000 )
double GiveAR(const double k){

return k*PixelWidth*iWidth/2000.0 ;

}

/* find such AR for internal LCM/J and LSM that level curves croses critical point zcr0 and it's preimages
for attracting ( also weakly attracting = parabolic) periodic point za0

it may fail if one iteration is bigger then smallest distance between periodic point and Julia set
*/
double GiveTunedAR(int i_Max){

complex double z= zcr0; // criical point
int i;
//int i_Max = 1000;
// critical point escapes very fast here. Higher valus gives infinity
for (i=0; i< i_Max; ++i ){
z= z*z*z +c*z; // forward iteration
}
double r = cabs(z-za0);
if ( r > AR_max ) {r = AR_max;}

return r;

}


Steps for parabolic basin

• choose component with critical point inside
• choose trap

Trap is a disk

• inside component with critical point inside
• trap has parabolic point on it's boundary
• center of the trap is the midpoint between last point of critical orbit and fixed point
• radius of the trap is half of distance between fixed point and last point of critical orbit

### Decomposition of target set

Examples

#### Binary decomposition

Here color of pixel ( exterior of Julia set) is proportional to sign of imaginary part of last iteration ( cimag) = the radial borders of are at dyadic angles (shows external rays at dyadic angles).

Main loop is the same as in escape time.

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

${\displaystyle T_{b}+=\{z:abs(z_{n})>ER~~{\mbox{and}}~~im(z_{n})>0\}\,}$

${\displaystyle 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;


Description

unsigned char ComputeColorOfBD (complex double z)
{

double cabsz;

int i;			// number of iteration
for (i = 0; i < IterMax_LSM; ++i)
{

cabsz = cabs(z); // numerical speed up : cabs(zp-z) = cabs(z) because zp = zcr = 0

//
if ( cabsz > ER  ||  cabsz < AR ) // if z is inside target set ( orbit trap)
{
if (cimag(z) > 0) // binary decomposition of target set
{  return 0;}
else {return 255; }

}

z = f(z);

}

return iColorOfUnknown;

}

##### rotation of target set by internal angle
 // for MBD
double t0 = 1.0 / 3.0; // period = 3

// Modified BD
unsigned char ComputeColorOfMBD (complex double z)
{

double cabsz;
double turn;

int i;			// number of iteration
for (i = 0; i < IterMax_LSM; ++i)
{

cabsz = cabs(z); // numerical speed up : cabs2(zp-z) = cabs2(z) because zp = zcr = 0

//  if z is inside target set ( orbit trap) = exterior of circle with radius ER
if ( cabsz > ER  ) // exterior
{
if (creal(z) > 0) // binary decomposition of target set
{  return 0;}
else {return 255; }

}

if ( cabsz  < AR ) // if z is inside target set ( orbit trap) = interior of circle with radius AR
{
turn = c_turn(z);
if (turn < t0 || turn > t0+0.5) // modified binary decomposition of target set
{  return 0;}
else {return 255; }

}

z = f(z);

}

return iColorOfUnknown;

}

##### Constant number of iterations

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.

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

## DEM/J

This algorithm has 2 versions:

Compare it with version for parameter plane and Mandelbrot set : DEM/M It's the same as M-set exterior distance estimation but with derivative w.r.t. Z instead of w.r.t. C.

## Convergence

In this algorithm distances between 2 points of the same orbit are checked

### average discrete velocity of orbit

average discrete velocity of orbit - code and description

It is used in case of :

### Cauchy Convergence Algorithm (CCA)

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

## Normality

Normality In this algorithm distances between points of 2 orbits are checked

### Checking equicontinuity by Michael Becker

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

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

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

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

## Koenigs coordinate

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

# Optimisation

## Distance

You don't need a square root to compare distances.[22]

## Symmetry

Julia sets can have many symmetries [23][24]

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 :[25]

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


Algorithm :

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

# Sets

## Target set

Target set or trap

One can divide it according to :

• attractors ( finite or infinite)
• dynamics ( hyperbolic, parabolic, elliptic )
• shape ( bailout test)
• decomposition of target set: binary decomposition ( BDM), angular decomposition,

### For infinite attractor - hyperbolic case

Target set ${\displaystyle 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 ${\displaystyle z=0\,}$ and radius =ER :

${\displaystyle 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 ${\displaystyle s\,}$ centered at origin

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

#### Julia sets

Escher like tilings is a modification of the level set method ( LSM/J). Here Level sets of escape time are different because targest set is different. Here target set is a scalled filled-in Julia set.

For more description see

• Fractint : escher_julia
• page 187 from The Science of fractal images by Heinz-Otto Peitgen, Dietmar Saupe, Springer [27]

### For finite attractors

internal level sets around fixed point

See :

## 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 )[30]

• 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 scanning [31]
• edge detection

## Fatou set

Interior of filled Julia set can be coloured :

More is here

# Video

One can make videos using :

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

Examples :

# References

1. Standard coloring algorithms from Ultra Fractal
2. new fractalforum : lowest-optimal-bailout-values-for-the-mandelbrot-sets/
4. Faster Fractals Through Algebra by Bruce Dawson ( author of Fractal eXtreme)
5. C code with gsl from tensorpudding
6. Euler examples by R. Grothmann
7. Drawing the Mandelbrot set by the method of escape lines. M. Romera et al.
8. Julia Curves, Mandelbrot Set, Harold V. McIntosh.
9. PythonDataScienceHandbook: density-and-contour-plots by Jake VanderPlas
10. math.stackexchange question: what-do-level-curves-signify
11. Contour lines by Rodolphe Vaillant
12. The fixed points and periodic orbits by E Demidov
13. Video : Julia Set Morphing with Magnetic Field lines by bryceguy72
14. Video : Mophing Julia set with color bands / stripes by FreymanArt
15. Gerard Westendorp : Platonic tilings of Riemann surfaces - 8 times iterated Automorphic function z->z^2 -0.1+ 0.75i
16. 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
17. Joseph H. Silverman : The arithmetic of dynamical systems. Springer, 2007. ISBN 0387699031, 9780387699035; page 22
18. Visualising Julia sets by Georg-Johann
19. Problem : How changes distance between 2 near points under iteration ? Can I tell to which set these points belong when I know it ?
20. Julia sets by Michael Becker. See the metric d(z,w)
21. Algorithms : Distance_approximations in wikibooks
22. The Julia sets symmetry by Evgeny Demidov
23. mathoverflow : symmetries-of-the-julia-sets-for-z2c
24. htJulia Jewels: An Exploration of Julia Sets by Michael McGoodwin (March 2000)
25. julia sets in Matlab by Jonas Lundgren
26. {Peitgen, H.O. and Fisher, Y. and Saupe, D. and McGuire, M. and Voss, R.F. and Barnsley, M.F. and Devaney, R.L. and Mandelbrot, B.B.}, (2012). The Science of Fractal Images. Springer Science & Business Media, 2012. p. 187. ISBN 9781461237846.
28. Tessellation of the Interior of Filled Julia Sets by Tomoki Kawahira
29. Mark Braverman : On efficient computation of parabolic Julia sets
30. ALGORITHM OF COMPUTER MODELLING OF JULIA SET IN CASE OF PARABOLIC FIXED POINT N.B.Ampilova, E.Petrenko
31. Ray Tracing Quaternion Julia Sets on the GPU by Keenan Crane
32. Tessellation of the Interior of Filled Julia Sets by Tomoki Kawahira
33. Julia-Set-Animations at devianart