Fractals/Iterations in the complex plane/Julia set

"... a single algorithm for computing all quadratic Julia sets does not exist."

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



=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) ?
 * escape radius

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

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

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

C with complex type from GSL :

C++ versions:

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

where :

Delphi version without explicit definition of complex numbers :

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

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.

Normalized iteration count (real escape time or fractional iteration or Smooth Iteration Count Algorithm (SICA))
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);

description:
 * FF : julia-smooth-colouring-how-to-do
 * stefan bion : fraktal-generator/colormapping/
 * FF smooth-histogram-rendering/
 * FF creating-a-good-palette-using-bezier-interpolation/

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

Density of level curves

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

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

How to control level curves
 * escape radius
 * shape of target set
 * manually
 * drawing equipotential lines
 * changing level sets ( level curves are boundaries of level sets)

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 ?



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

Internal Level Sets
See :
 * algorithm 0 of program Mandel by Wolf Jung

How to choose size of attracting trap such that level curves cross at critical point ?

It depends on
 * dynamic type ( superattracting/attracting, parabolic, repelling)
 * period ( child period in the parabolic case)


 * petal in the parabolic case :
 * for period 1 and 2 : radius of a circle with parabolic point on it's boundary
 * for higher periods sector of the circle with center at parabolic point
 * for other cases ( except repelling) it is radius of the circle with attractor as a center

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
 * fractint decomposition
 * fractalzoomer Binary_Decomposition

Binary decomposition
These curves :
 * are boundaries of binare decomposition boxes
 * are not field lines of potential = external rays

Notes
 * if the escape radius is too low, then binary (or ternary, etc) decomposition rays will have visible discontinuities at iteration bands. increasing escape radius makes discontinuities smaller, but changes the aspect ratio
 * mrob says exp(pi) is the best escape radius for binary decomposition as it makes boxes have square aspect ratio (possibly more visible when using exponential map transformation?)

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.

It creates radial level sets.

See also:
 * mandel : algorithm 9 =  zeros of qn(c)
 * video by bryceguy72
 * video by FreymanArt

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

BDM in the attracting domain
BDM in the attracting domain ( usually interior of Julia sets ) gives (pseudo) Field lines

Explanation by Gert Buschmann

In each Fatou domain (that is not neutral) there are two systems of lines orthogonal to each other: the equipotential lines (for the potential function or the real iteration number) and the field lines.

If we colour the Fatou domain according to the iteration number (and not the real iteration number $$\nu(z)$$, as defined in the previous section), the bands of iteration show the course of the equipotential lines. If the iteration is towards ∞ (as is the case with the outer Fatou domain for the usual iteration $$z^{2} + c$$), we can easily show the course of the field lines, namely by altering the colour according as the last point in the sequence of iteration is above or below the x-axis (first picture), but in this case (more precisely: when the Fatou domain is super-attracting) we cannot draw the field lines coherently - at least not by the method we describe here. In this case a field line is also called an external ray.

Let z be a point in the attracting Fatou domain. If we iterate z a large number of times, the terminus of the sequence of iteration is a finite cycle C, and the Fatou domain is (by definition) the set of points whose sequence of iteration converges towards C. The field lines issue from the points of C and from the (infinite number of) points that iterate into a point of C. And they end on the Julia set in points that are non-chaotic (that is, generating a finite cycle). Let r be the order of the cycle C (its number of points) and let z* be a point in C. We have $$f(f(\dots f(z^*))) = z^*$$ (the r-fold composition), and we define the complex number α by


 * $$\alpha = (d(f(f(\dots f(z))))/dz)_{z=z^*}.$$

If the points of C are $$z_i, i = 1, \dots, r (z_1 = z^*)$$, α is the product of the r numbers $$f'(z_i)$$. The real number 1/|α| is the attraction of the cycle, and our assumption that the cycle is neither neutral nor super-attracting, means that 1 < 1/|α| < ∞. The point z* is a fixed point for $$f(f(\dots f(z)))$$, and near this point the map $$f(f(\dots f(z)))$$ has (in connection with field lines) character of a rotation with the argument β of α (that is, $$\alpha = |\alpha|e^{\beta i}$$).

In order to colour the Fatou domain, we have chosen a small number ε and set the sequences of iteration $$z_k (k = 0, 1, 2, \dots, z_0 = z)$$ to stop when $$|z_k - z^*| < \epsilon$$, and we colour the point z according to the number k (or the real iteration number, if we prefer a smooth colouring). If we choose a direction from z* given by an angle θ, the field line issuing from z* in this direction consists of the points z such that the argument ψ of the number $$z_k - z^*$$ satisfies the condition that


 * $$\psi - k\beta = \theta \mod \pi. \, $$

For if we pass an iteration band in the direction of the field lines (and away from the cycle), the iteration number k is increased by 1 and the number ψ is increased by β, therefore the number $$\psi - k\beta \mod \pi$$ is constant along the field line.



A colouring of the field lines of the Fatou domain means that we colour the spaces between pairs of field lines: we choose a number of regularly situated directions issuing from z*, and in each of these directions we choose two directions around this direction. As it can happen that the two field lines of a pair do not end in the same point of the Julia set, our coloured field lines can ramify (endlessly) in their way towards the Julia set. We can colour on the basis of the distance to the center line of the field line, and we can mix this colouring with the usual colouring. Such pictures can be very decorative (second picture).

A coloured field line (the domain between two field lines) is divided up by the iteration bands, and such a part can be put into a one-to-one correspondence with the unit square: the one coordinate is (calculated from) the distance from one of the bounding field lines, the other is (calculated from) the distance from the inner of the bounding iteration bands (this number is the non-integral part of the real iteration number). Therefore, we can put pictures into the field lines (third picture).

ToDo

 * add slope to white

Inverse Iteration Method (IIM/J) : Julia set
Inverse iteration of repellor for drawing Julia set

Complex potential - Boettcher coordinate
See description here

DEM/J
This algorithm has 2 versions:
 * exterior
 * interior

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
It is used in case of :
 * parabolic dynamic
 * eliptic dynamic ( Siegel disc )

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)

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.

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.

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=
 * color in computer graphic
 * Color Theory/Color gradient
 * 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
 * Technicolor Julias by Shawn Hargreaves

=Sets=

Target set
Target set
 * trap for forward orbit
 * it is a set which captures any orbit tending to fixed / periodic point.

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 scanning
 * 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 converted to colour ( or shade of gray )
 * Internal Level Sets
 * attracting time ( sth like escape time but checks if (abs(z-attractor)<Attracting_radius
 * Binary decomposition
 * Tessellation of the Interior of Filled Julia Sets by Tomoki KAWAHIRA
 * discrete velocity in Siegel disc case

Periodic points
More is here

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

Examples :
 * Target set for internal ray 0 video
 * Quadratic Julia set with Internal level sets for internal ray 0 video

=More tutorials and code=
 * hypercomplex iterations - book
 * hvidtfeldts : distance-estimated-3d-fractals-v-the-mandelbulb-different-de-approximations/
 * 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
 * "Spiral Structures in Julia Sets and Related Sets", M. Michelitsch and O. E. Roessler in a book : SPIRAL SYMMETRY I. Hargittai and C. Pickover. (1992) World Scientific Publishing,
 * The Evolution of a Three-armed Spiral in the Julia Set, and Higher Order Spirals", A. G. Davis Philip in a book : SPIRAL SYMMETRY I. Hargittai and C. Pickover. (1992) World Scientific Publishing,
 * Beardon, A. : Symmetries of julia sets. The Mathematical Intelligencer. 1996-03-01 Springer New York Issn: 0343-6993 page 43 - 44.