Fractals/Iterations in the complex plane/qiterations
Iteration in mathematics refer to the process of iterating a function i.e. applying a function repeatedly, using the output from one iteration as the input to the next.^{[1]} Iteration of apparently simple functions can produce complex behaviours and difficult problems.^{[2]}
Applications
editOne can make inverse ( backward iteration) :
 of repeller for drawing Julia set ( IIM/J)^{[3]}
 of circle outside Jlia set (radius=ER) for drawing level curves of escape time ( which tend to Julia set)^{[4]}
 of circle inside Julia set (radius=AR) for drawing level curves of attract time ( which tend to Julia set)^{[5]}
 of critical orbit ( in Siegel disc case) for drawing Julia set ( probably only in case of Goldem Mean )
 for drawing external ray
Repellor for forward iteration is attractor for backward iteration
Notes
edit Iteration is always on the dynamic plane.
 There is no dynamic on the parameter plane.
 Mandelbrot set carries no dynamics. It is a set of parameter values.
 There are no orbits on parameter plane, one should not draw orbits on parameter plane.
 Orbit of critical point is on the dynamical plane
Iteration theory
editIt is a section from Tetration forum by Andrew Robbins 20060215 by Andrew Robbins
"Iteration is fundamental to dynamics, chaos, analysis, recursive functions, and number theory. In most cases the kind of iteration required in these subjects is integer iteration, i.e. where the iteration parameter is an integer. However, in the study of dynamical systems continuous iteration is paramount to the solution of some systems.
Different kinds of iteration can be classified as follows:
 Discrete Iteration
 Integer Iteration
 Fractional Iteration or Rational Iteration
 Nonanalytic Fractional Iteration
 analytic Fractional Iteration
 Continuous Iteration
Discrete iteration
editIterated function
 examples^{[6]}
 derivative^{[7]}
Integer iteration
editThe usual definition of iteration, where the functional equation:
is used to generate the sequence:
known as the natural iterates of f(x), which forms a monoid under composition.
For invertible functions f(x), the inverses are also considered iterates, and form the sequence:
known as the integer iterates of f(x), which forms a group under composition.
Fractional Iteration or Rational Iteration
editSolving the functional equation: . Once this functional equation is solved, then the rational iterates are the integer iterates of .
Nonanalytic Fractional Iteration
editBy choosing a nonanalytic fractional iterate, there is no uniqueness of the solutions obtained. (Iga's method)
Analytic Fractional Iteration
editBy solving for an analytic fractional iterate, there is a unique solution obtained in this way. (Dahl's method)
Continuous Iteration
editA generalization of the usual notion of iteration, where the functional equation (FE): f n(x) = f(f n1(x)) must be satisfied for all n in the domain (real or complex). If this is not the case, then to discuss these kinds of "iteration" (even though they are not technically "iteration" since they do not obey the FE of iteration), we will talk about them as though they were expressions for f n(x) where 0 ≤ Re(n) ≤ 1 and defined elsewhere by the FE of iteration. So even though a method is analytic, if it doesn't satisfy this fundamental FE, then by this redefinition, it is nonanalytic.
Nonanalytic Continuous Iteration
editBy choosing a nonanalytic continuous iterate, there is no uniqueness of the solutions obtained. (Galidakis' and Woon's methods)
Analytic Continuous Iteration or just Analytic Iteration
editBy solving for an analytic continuous iterate, there is a unique solution obtained in this way.
Realanalytic Iteration
editComplexanalytic Iteration or Holomorphic Iteration
editStep
edit Integer
 Fractional
 Continuous Iteration of Dynamical Maps.^{[8]}^{[9]} This continuous iteration can be discretized by the finite element method and then solved—in parallel—on a computer.
visualisation
edit folding ^{[10]}^{[11]}^{[12]}^{[13]}
 complex geometric sequence by Loo Kang Lawrence WEE
decomposition
editMove during iteration in case of complex quadratic polynomial is complex. It consists of 2 moves :
 angular move = rotation ( see doubling map)
 radial move ( see external and internal rays, invariant curves )
 fallin into target set and attractor ( in hyperbolic and parabolic case )

Radius abs(z) = r < 1.0 after some iterations using f0(z) = z*z

Distance between repetive iteration of point smaller then one

angle 1/7 after doubling map

angle 1/15 after doubling map

angle 11/36 after doubling map
angular move (rotation)
editCompute argument in turns^{[14]} of the complex number :
/*
gives argument of complex number in turns
*/
double GiveTurn( double complex z){
double t;
t = carg(z);
t /= 2*pi; // now in turns
if (t<0.0) t += 1.0; // map from (1/2,1/2] to [0, 1)
return (t);
}
direction
editforward
editbackward
editBackward iteration or inverse iteration^{[15]}
 Peitgen
 W Jung
 John Bonobo^{[16]}
Peitgen
edit /* Zn*Zn=Z(n+1)c */
Zx=ZxCx;
Zy=ZyCy;
/* sqrt of complex number algorithm from Peitgen, Jurgens, Saupe: Fractals for the classroom */
if (Zx>0)
{
NewZx=sqrt((Zx+sqrt(Zx*Zx+Zy*Zy))/2);
NewZy=Zy/(2*NewZx);
}
else /* ZX <= 0 */
{
if (Zx<0)
{
NewZy=sign(Zy)*sqrt((Zx+sqrt(Zx*Zx+Zy*Zy))/2);
NewZx=Zy/(2*NewZy);
}
else /* Zx=0 */
{
NewZx=sqrt(fabs(Zy)/2);
if (NewZx>0) NewZy=Zy/(2*NewZx);
else NewZy=0;
}
};
if (rand()<(RAND_MAX/2))
{
Zx=NewZx;
Zy=NewZy;
}
else {Zx=NewZx;
Zy=NewZy; }
Mandel
editHere is example of c code of inverse iteration using code from program Mandel by Wolf Jung
/*
/*
gcc i.c lm Wall
./a.out
z = 0.000000000000000 +0.000000000000000 i
z = 0.229955135116281 0.141357981605006 i
z = 0.378328716195789 0.041691618297441 i
z = 0.414752103217922 +0.051390827017207 i
*/
#include <stdio.h>
#include <math.h> // M_PI; needs lm also
#include <complex.h>
/* find c in component of Mandelbrot set
uses code by Wolf Jung from program Mandel
see function mndlbrot::bifurcate from mandelbrot.cpp
http://www.mndynamics.com/indexp.html
*/
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int Period)
{
//0 <= InternalRay<= 1
//0 <= InternalAngleInTurns <=1
double t = InternalAngleInTurns *2*M_PI; // from turns to radians
double R2 = InternalRadius * InternalRadius;
double Cx, Cy; /* C = Cx+Cy*i */
switch ( Period ) // of component
{
case 1: // main cardioid
Cx = (cos(t)*InternalRadius)/2(cos(2*t)*R2)/4;
Cy = (sin(t)*InternalRadius)/2(sin(2*t)*R2)/4;
break;
case 2: // only one component
Cx = InternalRadius * 0.25*cos(t)  1.0;
Cy = InternalRadius * 0.25*sin(t);
break;
// for each iPeriodChild there are 2^(iPeriodChild1) roots.
default: // higher periods : to do, use newton method
Cx = 0.0;
Cy = 0.0;
break; }
return Cx + Cy*I;
}
/* mndyncxmics::root from mndyncxmo.cpp by Wolf Jung (C) 20072014. */
// input = x,y
// output = u+v*I = sqrt(x+y*i)
complex double GiveRoot(complex double z)
{
double x = creal(z);
double y = cimag(z);
double u, v;
v = sqrt(x*x + y*y);
if (x > 0.0)
{ u = sqrt(0.5*(v + x)); v = 0.5*y/u; return u+v*I; }
if (x < 0.0)
{ v = sqrt(0.5*(v  x)); if (y < 0.0) v = v; u = 0.5*y/v; return u+v*I; }
if (y >= 0.0)
{ u = sqrt(0.5*y); v = u; return u+v*I; }
u = sqrt(0.5*y);
v = u;
return u+v*I;
}
// from mndlbrot.cpp by Wolf Jung (C) 20072014. part of Madel 5.12
// input : c, z , mode
// c = cx+cy*i where cx and cy are global variables defined in mndynamo.h
// z = x+y*i
//
// output : z = x+y*i
complex double InverseIteration(complex double z, complex double c, char key)
{
double x = creal(z);
double y = cimag(z);
double cx = creal(c);
double cy = cimag(c);
// f^{1}(z) = inverse with principal value
if (cx*cx + cy*cy < 1e20)
{
z = GiveRoot(x  cx + (y  cy)*I); // 2nd inverse function = key b
if (key == 'B') { x = x; y = y; } // 1st inverse function = key a
return z;
}
//f^{1}(z) = inverse with argument adjusted
double u, v;
complex double uv ;
double w = cx*cx + cy*cy;
uv = GiveRoot(cx/w (cy/w)*I);
u = creal(uv);
v = cimag(uv);
//
z = GiveRoot(w  cx*x  cy*y + (cy*x  cx*y)*I);
x = creal(z);
y = cimag(z);
//
w = u*x  v*y;
y = u*y + v*x;
x = w;
if (key =='A'){
x = x;
y = y; // 1st inverse function = key a
}
return x+y*I; // key b = 2nd inverse function
}
/*f^{1}(z) = inverse with argument adjusted
"When you write the real and imaginary parts in the formulas as complex numbers again,
you see that it is sqrt( c / c^2 ) * sqrt( c^2  conj(c)*z ) ,
so this is just sqrt( z  c ) except for the overall sign:
the standard squareroot takes values in the right halfplane, but this is rotated by the squareroot of c .
The new line between the two planes has half the argument of c .
(It is not orthogonal to c ... )"
...
"the argument adjusting in the inverse branch has nothing to do with computing external arguments. It is related to itineraries and kneading sequences, ...
Kneading sequences are explained in demo 4 or 5, in my slides on the stripping algorithm, and in several papers by Bruin and Schleicher.
W Jung " */
double complex GiveInverseAdjusted (complex double z, complex double c, char key){
double t = cabs(c);
t = t*t;
z = csqrt(c/t)*csqrt(tz*conj(c));
if (key =='A') z = z; // 1st inverse function = key a
// else key == 'B'
return z;
}
// make iMax inverse iteration with negative sign ( a in Wolf Jung notation )
complex double GivePrecriticalA(complex double z, complex double c, int iMax)
{
complex double za = z;
int i;
for(i=0;i<iMax ;++i){
printf("i = %d , z = (%f, %f) \n ", i, creal(za), cimag(za) );
za = InverseIteration(za,c, 'A');
}
printf("i = %d , z = (%f, %f) \n ", i, creal(za), cimag(za) );
return za;
}
// make iMax inverse iteration with negative sign ( a in Wolf Jung notation )
complex double GivePrecriticalA2(complex double z, complex double c, int iMax)
{
complex double za = z;
int i;
for(i=0;i<iMax ;++i){
printf("i = %d , z = (%f, %f) \n ", i, creal(za), cimag(za) );
za = GiveInverseAdjusted(za,c, 'A');
}
printf("i = %d , z = (%f, %f) \n ", i, creal(za), cimag(za) );
return za;
}
int main(){
complex double c;
complex double z;
complex double zcr = 0.0; // critical point
int iMax = 10;
int iPeriodChild = 3; // period of
int iPeriodParent = 1;
c = GiveC(1.0/((double) iPeriodChild), 1.0, iPeriodParent); // root point = The unique point on the boundary of a muatom of period P where two external angles of denominator = (2^P1) meet.
z = GivePrecriticalA( zcr, c, iMax);
printf("iAngle = %d/%d c = (%f, %f); z = (%.16f, %.16f) \n ", iPeriodParent, iPeriodChild, creal(c), cimag(c), creal(z), cimag(z) );
z = GivePrecriticalA2( zcr, c, iMax);
printf("iAngle = %d/%d c = (%f, %f); z = (%.16f, %.16f) \n ", iPeriodParent, iPeriodChild, creal(c), cimag(c), creal(z), cimag(z) );
return 0;
}
Test
editOne can iterate ad infinitum. Test tells when one can stop
 bailout test for forward iteration
Target set or trap
editTarget set is used in test. When zn is inside target set then one can stop the iterations.
Planes
editParameter plane
edit"Mandelbrot set carries no dynamics. It is a set of parameter values. There are no orbits on parameter plane, one should not draw orbits on parameter plane. Orbit of critical point is on the dynamical plane"
Dynamic plane
edit"The polynomial Pc maps each dynamical ray to another ray doubling the angle (which we measure in full turns, i.e. 0 = 1 = 2π rad = 360◦), and the dynamical rays of any polynomial “look like straight rays” near infinity. This allows us to study the Mandelbrot and Julia sets combinatorially, replacing the dynamical plane by the unit circle, rays by angles, and the quadratic polynomial by the doubling modulo one map." Virpi K a u k o^{[17]}
Dynamic plane f0 for c=0
editLets take c=0, then one can call dynamical plane plane.
Here dynamical plane can be divided into :
 Julia set =
 Fatou set which consists of 2 subsets :
 interior of Julia set = basin of attraction of finite attractor =
 exterior of Julia set = basin of attraction of infinity =
Forward iteration
edit
where :
 r is the absolute value or modulus or magnitude of a complex number z = x + i
 is the argument of complex number z (in many applications referred to as the "phase") is the angle of the radius with the positive real axis. Usually principal value is used
so
and forward iteration :^{[18]}
Forward iteration:
 squares radius and doubles angle ( phase, argument)^{[19]}^{[20]}
 gives forward orbit = list of points {z0, z1, z2, z3... , zn} which lays on exponential spirals.^{[21]}^{[22]}
One can check it interactively :
 Geogebra plet by Jorj KowszunLecturer  try move along radial lines from the center to see changing behaviour
 Visualization of complex_square_iteration by M McClure
 program Mandel by W Jung
 squaring complex numbers by Jeffrey Ventrella
 MIT mathlets : complexexponential
Chaos and the complex squaring map
editThe informal reason why the iteration is chaotic is that the angle doubles on every iteration and doubling grows very quickly as the angle becomes ever larger, but angles which differ by multiples of 2π radians are identical. Thus, when the angle exceeds 2π, it must wrap to the remainder on division by 2π. Therefore, the angle is transformed according to the dyadic transformation (also known as the 2x mod 1 map). As the initial value z_{0} has been chosen so that its argument is not a rational multiple of π, the forward orbit of z_{n} cannot repeat itself and become periodic.
More formally, the iteration can be written as:
where is the resulting sequence of complex numbers obtained by iterating the steps above, and represents the initial starting number. We can solve this iteration exactly:
Starting with angle θ, we can write the initial term as so that . This makes the successive doubling of the angle clear. (This is equivalent to the relation .)
Escape test
editIf distance between:
 point z of exterior of Julia set
 Julia set
is :
then point z escapes (= it's magnitude is greate then escape radius = ER):
after :
 steps in nonparabolic case
 steps in parabolic case ^{[23]}
See also:
Backward iteration
editEvery angle (real number) α ∈ R/Z measured in turns has :
 one image = 2α mod 1 under doubling map
 "two preimages under the doubling map: ".^{[24]} Inverse of doubling map is multivalued function.
Note that difference between these 2 preimages
is half a turn = 180 degrees = Pi radians.
On complex dynamical plane backward iteration using quadratic polynomial
gives backward orbit = binary tree of preimages :
One can't choose good path in such tree without extra informations.
Not that preimages show rotational symmetry ( 180 degrees)
For other functions see Fractalforum^{[25]}
See also:
Dynamic plane for fc
editOne can check it with :
 M McClure : visualization : interactive_basins
 GeoGebra  Complex Iteration  The Mandelbrot Set by Mark Dabbs
 critical orbits by Jeremy Rifkin
Level curves of escape or attracting time
edit
Preimages of circle from the exterior gives level sets of the escaping time

Preimages of circle from the interior gives level sets of the attracting time
Julia set by IIM/J
editTheory
 the periodic points are dense in the Julia set
 Julia set is the closure of the set of repelling periodic points
So drawing repelling periodic point and it's orbit ( forward and backward= inverse) gives visually good aproximation of Julia set = set of points dense enough that nonuniform distribution of these points over Julia set is not important.

Preimages of repelling fixed point tend to Julia set for C = i. Image and source code

Julia set made with MIIM. Image and maxima source code

Preimages of critical orbit tend to whole Julia set in Siegel disc case

distribution of points in simple IIM/J

c = 0.750357820200574 +0.047756163825227 i

C = ( 0.4 0.3 )
In escape time one computes forward iteration of point z.
In IIM/J one computes:
 repelling fixed point^{[26]} of complex quadratic polynomial
 preimages of by inverse iterations
"We know the periodic points are dense in the Julia set, but in the case of weird ones (like the ones with Cremer points, or even some with Siegel disks where the disk itself is very 'deep' within the Julia set, as measured by the external rays), the periodic points tend to avoid certain parts of the Julia set as long as possible. This is what causes the 'inverse method' of rendering images of Julia sets to be so bad for those cases." Jacques Carette^{[27]}
Because square root is multivalued function then each has two preimages . Thus inverse iteration creates binary tree of preimages.

binary tree

expanded growth of binary tree
Because of expanded growth of binary tree, the number of preimages grows exponentialy: the number of nodes in a full binary tree is
where
 is the height of the tree
If one wants to draw full binary tree then the methods of storing binary trees can waste a fair bit of memory, so alternative is
 threaded binary tree
 draw only some path from binary tree,
 random path : the longest path = path from roottoleaf
See also :
 online tool : Wolfram Alpha^{[28]}^{[29]}
Root of tree
edit repelling fixed point^{[30]} of complex quadratic polynomial
  beta
 other repelling periodic points ( cut points of filled Julia set ). It will be important especially in case of the parabolic Julia set.
"... preimages of the repelling fixed point beta. These form a tree like
beta beta beta beta beta x y
So every point is computed at last twice when you start the tree with beta. If you start with beta, you will get the same points with half the number of computations.
Something similar applies to the preimages of the critical orbit. If z is in the critical orbit, one of its two preimages will be there as well, so you should draw z and the tree of its preimages to avoid getting the same points twice." (Wolf Jung )
Variants of IIM
edit random choose one of two roots IIM ( up to chosen level max). Random walk through the tree. Simplest to code and fast, but inefficient. Start from it.
 both roots with the same probability
 more often one then other root
 draw all roots ( up to chosen level max)^{[31]}
 using recurrence
 using stack ( faster ?)
 draw some rare paths in binary tree = MIIM. This is modification of drawing all roots. Stop using some rare paths.
 using hits table ( while hit(pixel(iX,iY)) < HitMax )^{[32]}
 using derivative ( while derivative(z) < DerivativeMax)^{[33]}^{[34]}:
 " speed up the calculation by pruning dense branches of the tree. One such method is to prune branches when the map becomes contractive (the cumulative derivative becomes large). " Paul Nylander^{[35]}
 " We cut off the sub tree from a given if the derivative is greater then limit. This eliminates dominant highly contractive regions of the inverse iteration, which have already been registered. We can calculate successive derivatives iteratively. Colour by the log of the absolute derivative" ^{[36]}
Examples of code :
Compare it with:
See also
 color: parts of Julia set have different color
See also
edit Dynamical systems
 Fixed points
 Lyapunov number
 Functional equations
 Abel function
 Schroeder function
 Boettcher function
 Julia function
 Special matrices
 Carleman matrix
 Bell matrix
 AbelRobbins matrix
References
edit ↑ wikipedia : Iteration
 ↑ From local to global theories of iteration by Vaclav Kucera
 ↑ Inverse Iteration Algorithms for Julia Sets by Mark McClure
 ↑ Complex iteration by Microcomputadoras
 ↑ On rational maps with two critical points by John Milnor, fig. 5
 ↑ Iterated Functions by Tom Davis
 ↑ The Existence and Uniqueness of the Taylor Series of Iterated Functions by Daniel Geisler
 ↑ Continuous Iteration of Dynamical Maps R. Aldrovandi, L. P. Freitas (Sao Paulo, IFT)
 ↑ Continuous_iteration_of_fractals by Gerard Westendorp
 ↑ howtofoldajuliafractal by Steven Wittens
 ↑ Folding a Circle into a Julia Set by Karl Sims
 ↑ Visual Explanation of the Complexity in Julia Sets by Okke Schrijvers, Jarke J. van Wijk ( see video in the supporting info)
 ↑ How to Build a Julia Set
 ↑ turn
 ↑ Understanding Julia and Mandelbrot Sets by Karl Sims
 ↑ Une méthode rapide pour tracer les ensembles de Julia : l'itération inverse by John Bonobo
 ↑ Trees of visible components in the Mandelbrot set by Virpi K a u k o , FUNDAMENTA MATHEMATICAE 164 (2000)
 ↑ Real and imaginary parts of polynomial iterates by Julia A. Barnes, Clinton P. Curry and Lisbeth E. Schaubroeck. New York J. Math. 16 (2010) 749–761.
 ↑ mandelbrothue by Richard A Thomson
 ↑ phase by Linas Vepstas
 ↑ Complex numbers by David E Joyce
 ↑ Powers of complex numbers from Suitcase of Dreams
 ↑ Parabolic Julia Sets are Polynomial Time Computable by Mark Braverman
 ↑ SYMBOLIC DYNAMICS AND SELFSIMILAR GROUPS by VOLODYMYR NEKRASHEVYCH
 ↑ Query about general Julia set IFS for higher powers.
 ↑ wikipedia : repelling fixed point
 ↑ mathbbc/185430 mathoverflow question: clusteringofperiodicpointsforapolynomialiterationofmathbbc
 ↑ Wolfram Alpha
 ↑ example
 ↑ wikipedia : repelling fixed point
 ↑ Fractint documentation  Inverse Julias
 ↑ Image and c source code : IIMM using hit limit
 ↑ Exploding the Dark Heart of Chaos by Chris King
 ↑ Drakopoulos V., Comparing rendering methods for Julia sets, Journal of WSCG 10 (2002), 155161
 ↑ bugman123: Fractals
 ↑ dhushara : DarkHeart