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]

One 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

• 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

It is a section fromTetration forum by Andrew Robbins 2006-02-15 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.

Differentkinds of iteration can be classified as follows:

• Discrete Iteration
  • Integer Iteration
  • Fractional Iteration or Rational Iteration
    • Non-analytic Fractional Iteration
    • analytic Fractional Iteration
• Continuous Iteration

Iterated function

• examples^[6]
• derivative^[7]

The usual definition of iteration, where the functional equation:

${\displaystyle f^n(x)=f(f^{n-1}(x))}$

is used to generate the sequence:

${\displaystyle \{f(x),f^2(x),f^3(x),...\}}$

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:

${\displaystyle \{...,f^{-2}(x),f^{-1}(x),x,f(x),f^2(x),...\}}$

known as the integer iterates of f(x), which forms a group under composition.

Solving the functional equation:${\displaystyle f(x)=g^n(x)}$. Once this functional equation is solved, then the rational iterates${\displaystyle f^{(m/n)}(x)}$ are the integer iterates of${\displaystyle g(x)}$.

By choosing a non-analytic fractional iterate, there is no uniqueness of the solutions obtained. (Iga's method)

By solving for an analytic fractional iterate, there is a unique solution obtained in this way. (Dahl's method)

A generalization of the usual notion of iteration, where the functional equation (FE): f n(x) = f(f n-1(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 re-definition, it is non-analytic.

By choosing a non-analytic continuous iterate, there is no uniqueness of the solutions obtained. (Galidakis' and Woon's methods)

By solving for an analytic continuous iterate, there is a unique solution obtained in this way.

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

• folding^{[10][11][12][13]}
• complex geometric sequence by Loo Kang Lawrence WEE

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

• radial move
• 
  Radius abs(z) = r < 1.0 after some iterations using f0(z) = z*z
• 
  Distance between repetive iteration of point smaller then one

• angular move
• 
  angle 1/7 after doubling map
• 
  angle 1/15 after doubling map
• 
  angle 11/36 after doubling map

Compute argument in turns^[14] of the complex number :

/*gives argument of complex number in turns*/doubleGiveTurn(doublecomplexz){doublet;t=carg(z);t/=2*pi;// now in turnsif(t<0.0)t+=1.0;// map from (-1/2,1/2] to [0, 1)return(t);}

Backward iteration or inverse iteration^[15]

• Peitgen
• W Jung
• John Bonobo^[16]

/* Zn*Zn=Z(n+1)-c */Zx=Zx-Cx;Zy=Zy-Cy;/* 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);elseNewZy=0;}};if(rand()<(RAND_MAX/2)){Zx=NewZx;Zy=NewZy;}else{Zx=-NewZx;Zy=-NewZy;}

Here is example of c code of inverse iteration using code from programMandel by Wolf Jung

/*/* gcc i.c -lm -Wall ./a.outz = 0.000000000000000  +0.000000000000000 iz = -0.229955135116281  -0.141357981605006 iz = -0.378328716195789  -0.041691618297441 iz = -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*/doublecomplexGiveC(doubleInternalAngleInTurns,doubleInternalRadius,unsignedintPeriod){//0 <= InternalRay<= 1//0 <= InternalAngleInTurns <=1doublet=InternalAngleInTurns*2*M_PI;// from turns to radiansdoubleR2=InternalRadius*InternalRadius;doubleCx,Cy;/* C = Cx+Cy*i */switch(Period)// of component{case1:// main cardioidCx=(cos(t)*InternalRadius)/2-(cos(2*t)*R2)/4;Cy=(sin(t)*InternalRadius)/2-(sin(2*t)*R2)/4;break;case2:// only one componentCx=InternalRadius*0.25*cos(t)-1.0;Cy=InternalRadius*0.25*sin(t);break;// for each iPeriodChild  there are 2^(iPeriodChild-1) roots.default:// higher periods : to do, use newton methodCx=0.0;Cy=0.0;break;}returnCx+Cy*I;}/* mndyncxmics::root from mndyncxmo.cpp  by Wolf Jung (C) 2007-2014. */// input = x,y// output = u+v*I = sqrt(x+y*i)complexdoubleGiveRoot(complexdoublez){doublex=creal(z);doubley=cimag(z);doubleu,v;v=sqrt(x*x+y*y);if(x>0.0){u=sqrt(0.5*(v+x));v=0.5*y/u;returnu+v*I;}if(x<0.0){v=sqrt(0.5*(v-x));if(y<0.0)v=-v;u=0.5*y/v;returnu+v*I;}if(y>=0.0){u=sqrt(0.5*y);v=u;returnu+v*I;}u=sqrt(-0.5*y);v=-u;returnu+v*I;}// from mndlbrot.cpp  by Wolf Jung (C) 2007-2014. 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*icomplexdoubleInverseIteration(complexdoublez,complexdoublec,charkey){doublex=creal(z);doubley=cimag(z);doublecx=creal(c);doublecy=cimag(c);// f^{-1}(z) = inverse with principal valueif(cx*cx+cy*cy<1e-20){z=GiveRoot(x-cx+(y-cy)*I);// 2-nd inverse function = key bif(key=='B'){x=-x;y=-y;}// 1-st inverse function = key areturn-z;}//f^{-1}(z) =  inverse with argument adjusteddoubleu,v;complexdoubleuv;doublew=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;// 1-st inverse function = key a}returnx+y*I;// key b =  2-nd 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 square-root 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 " */doublecomplexGiveInverseAdjusted(complexdoublez,complexdoublec,charkey){doublet=cabs(c);t=t*t;z=csqrt(-c/t)*csqrt(t-z*conj(c));if(key=='A')z=-z;// 1-st inverse function = key a// else key == 'B'returnz;}// make iMax inverse iteration with negative sign ( a in Wolf Jung notation )complexdoubleGivePrecriticalA(complexdoublez,complexdoublec,intiMax){complexdoubleza=z;inti;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));returnza;}// make iMax inverse iteration with negative sign ( a in Wolf Jung notation )complexdoubleGivePrecriticalA2(complexdoublez,complexdoublec,intiMax){complexdoubleza=z;inti;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));returnza;}intmain(){complexdoublec;complexdoublez;complexdoublezcr=0.0;// critical pointintiMax=10;intiPeriodChild=3;// period ofintiPeriodParent=1;c=GiveC(1.0/((double)iPeriodChild),1.0,iPeriodParent);// root point = The unique point on the boundary of a mu-atom of period P where two external angles of denominator = (2^P-1) 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));return0;}

One can iterate ad infinitum. Test tells when one can stop

• bailout test for forward iteration

Target set is used in test. When zn is inside target set then one can stop the iterations.

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

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

Lets take c=0, then one can call dynamical plane${\displaystyle f_0}$ plane.

Here dynamical plane can be divided into :

• Julia set =${\displaystyle \{z:|z|=1\}}$
• 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 =${\displaystyle \{z:|z|>1\}}$

${\displaystyle z=r{e}^{i\theta }}$

where :

• r is theabsolute value ormodulus ormagnitude of a complex number z = x + i
• ${\displaystyle \theta }$ 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

${\displaystyle r=|z|=\sqrt{x^2+y^2}.}$

so

${\displaystyle f_0(z)=z^2=(r{e}^{i\theta }{)}^2=r^2{e}^{i2\theta }}$

and forward iteration :^[18]

${\displaystyle f_0^n(z)=r^{2^n}{e}^{i{2}^n\theta }}$

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 Kowszun-Lecturer - 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 : complex-exponential

The 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 mustwrap 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 valuez₀ has been chosen so that its argument is not a rational multiple of π, the forward orbit ofz_n cannot repeat itself and become periodic.

More formally, the iteration can be written as:

${\displaystyle z_{n+1}=z_n^2}$

where${\displaystyle z_n}$ is the resulting sequence of complex numbers obtained by iterating the steps above, and${\displaystyle z_0}$ represents the initial starting number. We can solve this iteration exactly:

${\displaystyle z_n=z_0^{2^n}}$

Starting with angle θ, we can write the initial term as${\displaystyle z_0=\exp (i\theta )}$ so that${\displaystyle z_n=\exp (i{2}^n\theta )}$. This makes the successive doubling of the angle clear. (This is equivalent to the relation${\displaystyle z_n=\cos (2^n\theta )+i\sin (2^n\theta )}$.)

If distance between:

• point z of exterior of Julia set
• Julia set${\displaystyle J_c}$

is :

${\displaystyle distance(z,J_c)=2^{-n}}$

then point z escapes (= it's magnitude is greate thenescape radius = ER):

${\displaystyle |z_n|>ER}$

after :

• ${\displaystyle n}$ steps in non-parabolic case
• ${\displaystyle 2^n}$ steps inparabolic case^[23]

See also:

• the precision needed for escape test
• bailout test
• escape time
• david bau : complex function viewer

Every angle (real number) α ∈ R/Z measured in turns has :

• one image = 2α mod 1 underdoubling map
• "two preimages underthe doubling map:${\displaystyle \frac{\alpha }{2},\frac{\alpha +1}{2}}$".^[24] Inverse of doubling map is multivalued function.

Note that difference between these 2 preimages

${\displaystyle \frac{\alpha }{2}-\frac{\alpha +1}{2}=\frac{1}{2}}$

is half a turn = 180 degrees = Pi radians.

Images and preimages under doubling map d

${\displaystyle \alpha }$	${\displaystyle d^1(\alpha )}$	${\displaystyle d^{-1}(\alpha )}$
${\displaystyle \frac{1}{2}}$	${\displaystyle \frac{1}{1}}$	${\displaystyle \left\{\frac{1}{4},\frac{3}{4}\right\}}$
${\displaystyle \frac{1}{3}}$	${\displaystyle \frac{2}{3}}$	${\displaystyle \left\{\frac{1}{6},\frac{4}{6}\right\}}$
${\displaystyle \frac{1}{4}}$	${\displaystyle \frac{1}{2}}$	${\displaystyle \left\{\frac{1}{8},\frac{5}{8}\right\}}$
${\displaystyle \frac{1}{5}}$	${\displaystyle \frac{2}{5}}$	${\displaystyle \left\{\frac{1}{10},\frac{6}{10}\right\}}$
${\displaystyle \frac{1}{6}}$	${\displaystyle \frac{1}{3}}$	${\displaystyle \left\{\frac{1}{12},\frac{7}{12}\right\}}$
${\displaystyle \frac{1}{7}}$	${\displaystyle \frac{2}{7}}$	${\displaystyle \left\{\frac{1}{14},\frac{4}{7}\right\}}$

On complex dynamical plane backward iteration using quadratic polynomial${\displaystyle f_c}$

${\displaystyle f_c(z)=z^2+c}$

gives backward orbit =binary tree of preimages :

${\displaystyle z}$

${\displaystyle -\sqrt{z-c},+\sqrt{z-c}}$

${\displaystyle -\sqrt{-\sqrt{z-c}-c},+\sqrt{-\sqrt{z-c}-c},-\sqrt{+\sqrt{z-c}-c},+\sqrt{+\sqrt{z-c}-c}}$

One can't choose good path in such tree withoutextra informations.

Not that preimages showrotational symmetry ( 180 degrees)

For other functions see Fractalforum^[25]

See also:

• Julia Set Grid by Steve Phelps. The mapping z=sqrt(z-c) on geogebra

One can check it with :

• M McClure : visualization : interactive_basins
• GeoGebra - Complex Iteration - The Mandelbrot Set by Mark Dabbs
• critical orbits by Jeremy Rifkin

• 
  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

Theory

• 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${\displaystyle Z_0={\beta }_c}$
• preimages of${\displaystyle Z_0}$ byinverse iterations

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

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

Becausesquare root ismultivalued function then each${\displaystyle Z_n}$ has two preimages${\displaystyle Z_{n-1}}$. Thus inverse iteration createsbinary tree of preimages.

• 
  binary tree
• 
  expanded growth of binary tree

Because of expanded growth of binary tree, the number of preimagesgrows exponentialy: the number of nodes${\displaystyle n}$ in a full binary tree is

${\displaystyle n=2^{h+1}-1}$

where

• ${\displaystyle h}$ 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 root-to-leaf

See also :

• online tool : Wolfram Alpha^[28][29]

• repelling fixed point^[30] of complex quadratic polynomial${\displaystyle Z_0={\beta }_c}$
• - beta
• other repelling periodic points ( cut points of filled Julia set ). It will be important especially in case ofthe 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 numberof computations.

Something similar applies to the preimages of the critical orbit. If z isin the critical orbit, one of its two preimages will be there as well, so youshould draw -z and the tree of its preimages to avoid getting the samepoints twice." (Wolf Jung )

• 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 than the 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 :

  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 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 /* random choose of one of two roots  */   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.
• Matlab code

 /* Maxima CAS code */ /* Gives points of backward orbit of z=repellor       */ GiveBackwardOrbit(c,repellor,zxMin,zxMax,zyMin,zyMax,iXmax,iYmax):=  block(   hit_limit:4, /* proportional to number of details and time of drawing */   PixelWidth:(zxMax-zxMin)/iXmax,   PixelHeight:(zyMax-zyMin)/iYmax,   /* 2D array of hits pixels . Hit > 0 means that point was in orbit */   array(Hits,fixnum,iXmax,iYmax), /* no hits for beginning */  /* choose repeller z=repellor as a starting point */  stack:[repellor], /*save repellor in stack */  /* save first point to list of pixels  */   x_y:[repellor],  /* reversed iteration of repellor */  loop,  /* pop = take one point from the stack */  z:last(stack),  stack:delete(z,stack),  /*inverse iteration - first preimage (root) */  z:finverseplus(z,c),  /* translate from world to screen coordinate */  iX:fix((realpart(z)-zxMin)/PixelWidth),  iY:fix((imagpart(z)-zyMin)/PixelHeight),  hit:Hits[iX,iY],  if hit<hit_limit      then     (    Hits[iX,iY]:hit+1,    stack:endcons(z,stack), /* push = add z at the end of list stack */    if hit=0 then x_y:endcons( z,x_y)    ),  /*inverse iteration - second preimage (root) */  z:-z, /* translate from world to screen coordinate, coversion to integer */  iX:fix((realpart(z)-zxMin)/PixelWidth),  iY:fix((imagpart(z)-zyMin)/PixelHeight),  hit:Hits[iX,iY],  if hit<hit_limit      then     (     Hits[iX,iY]:hit+1,     stack:endcons(z,stack), /* push = add z at the end of list stack to continue iteration */     if hit=0 then x_y:endcons( z,x_y)    ),   if is(not emptyp(stack)) then go(loop),  return(x_y) /* list of pixels in the form [z1,z2] */ )$

#octavem-file#IIMusingHittable#saveasaiimm_hit.minoctaveworkingdir#runoctaveandiimm_hit##stack-likeoperation#http://www.gnu.org/software/octave/doc/interpreter/Miscellaneous-Techniques.html#Miscellaneous-Techniques#octavecanresizearray#a=[];#while(condition)#a(end+1)=value;#"push"operation#a(end)=[];#"pop"operation##--------------generaloctavesettings----------clearall;moreoff;pkgloadimage;#imwritepkgloadmiscellaneous;#waitbar#---------------constandvar-----------------------------#definesomeglobalvarATEACHLEVELwhereyouwanttouseit(PascalCdeMills)#?forglobalvariablesonecan'tdefineandgiveinitialvalueatthesametime#integer(screen)coordinateiSide=1000ixMax=iSideiyMax=iSide#globalHitLimit;HitLimit=1#proportionaltonumberofdetailsandtimeofdrawingglobalHits;#2Darrayofpixels.Hit>0meansthatpointwasinorbitHits=zeros(iyMax,ixMax);#imageasa2Dmatrixof24bitcolorscodedfrom0.0to1.0globalMyImage;MyImage=zeros(iyMax,ixMax,3);#matrixfilledwith0.0(blackimage)=[rows,columns,color_depth]#world(float)coordinate-dynamical(Z)planeglobaldSide;globalZxMin;globalZxMax;globalZyMin;globalZyMax;globalz;globalPixelHeight;globalPixelWidth;#addvaluestoglobalvariablesdSide=1.25ZxMin=-dSideZxMax=dSideZyMin=-dSideZyMax=dSidePixelHeight=(ZyMax-ZyMin)/(iyMax-1)PixelWidth=(ZxMax-ZxMin)/(ixMax-1)#pseudostack=resizablearrayglobalStack;Stack=[];globalStackIndex;StackIndex=0;c=complex(.27327401,0.00756218);globalColor24White;Color24White=[1.0,1.0,1.0];#----------------functions------------------function[iY, iX] = f(z)#definesomeglobalvarATEACHLEVELwhereyouwanttouseit(PascalCdeMills)globalZxMin;globalZyMax;globalPixelWidth;globalPixelHeight;#translatefromworldtoscreencoordinateiX=int32((real(z)-ZxMin)/PixelWidth);iY=int32((ZyMax-imag(z))/PixelHeight);#invertyaxisendfunction;#plotpointwithintegercoorfinatetoarrayMyImagefunctionPlot( iY, iX , color)globalMyImage;MyImage(iY,iX,1)=color(1);MyImage(iY,iX,2)=color(2);MyImage(iY,iX,3)=color(3);endfunction;#push=putpointzonthestackfunctionpush(z)globalStack;globalStackIndex;StackIndex=size(Stack,2);#updateglobalvarStackIndex=StackIndex+1;#updateglobalvarStack(StackIndex)=z;#"push"zonpseudostackendfunction;#pop=takepointzfromthestackfunctionz=pop()globalStack;globalStackIndex;StackIndex=size(Stack,2);#updateglobalvarz=Stack(StackIndex);#popzfrompseudostackStack(StackIndex)=[];#makepseudostackshorterStackIndex=StackIndex-1;#updateglobalvarendfunction;functionCheckPoint(z)#errorishereglobalHits;globalHitLimit;globalColor24White;[iY,iX]=f(z);hit=Hits(iY,iX);if(hit<HitLimit)push(z);#(putzonthestack)tocontinueiterationif(hit<1)Plot(iY,iX,Color24White);endif;#plotHits(iY,iX)=hit+1;#updateHitstableendif;endfunction;#CheckPoint#--------------------main---------------------------------------#Determinetheunstablefixedpointbeta#http://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappingsbeta=0.5+sqrt(0.25-c)z=-betaCheckPoint(z);while(StackIndex>0)#ifstackisnotemptyz=pop();#takepointzfromthestack#computeits2preimageszand-z(inversefunctionofcomplexquadraticpolynomial)z=sqrt(z-c);#checkpoints:draw,putonthestacktocontinueiterationsCheckPoint(z);CheckPoint(-z);endwhile;#-------------------image--------------------------------------imread("a7.png");#firstloadanyexistingimagetoresolvebug:paniccrashusingimwrite:octave:magick/semaphore.c:525firstloadanyimageimage(MyImage);#DisplayImagename=strcat("iim",int2str(HitLimit)," .png");imwrite(MyImage,name);#saveimagetofile.thisrequirestheImageMagick"convert"utility.

See also

• color: parts of Julia set have different color
  • virtual math museum: julia set
  • fractalforums.org: Julia set

• Dynamical systems
  • Fixed points
  • Lyapunov number
• Functional equations
  • Abel function
  • Schroeder function
  • Boettcher function
  • Julia function
• Special matrices
  • Carleman matrix
  • Bell matrix
  • Abel-Robbins matrix

1. ↑wikipedia : Iteration
2. ↑From local to global theories of iteration by Vaclav Kucera
3. ↑Inverse Iteration Algorithms for Julia Sets by Mark McClure
4. ↑Complex iteration by Microcomputadoras
5. ↑On rational maps with two critical points by John Milnor, fig. 5
6. ↑Iterated Functions by Tom Davis
7. ↑The Existence and Uniqueness of the Taylor Series of Iterated Functions by Daniel Geisler
8. ↑Continuous Iteration of Dynamical Maps R. Aldrovandi, L. P. Freitas (Sao Paulo, IFT)
9. ↑Continuous_iteration_of_fractals by Gerard Westendorp
10. ↑how-to-fold-a-julia-fractal by Steven Wittens
11. ↑Folding a Circle into a Julia Set by Karl Sims
12. ↑Visual Explanation of the Complexity in Julia Sets by Okke Schrijvers, Jarke J. van Wijk ( see video in the supporting info)
13. ↑How to Build a Julia Set
14. ↑turn
15. ↑Understanding Julia and Mandelbrot Sets by Karl Sims
16. ↑Une méthode rapide pour tracer les ensembles de Julia : l'itération inverse by John Bonobo
17. ↑Trees of visible components in the Mandelbrot set by Virpi K a u k o , FUNDAMENTA MATHEMATICAE 164 (2000)
18. ↑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.
19. ↑mandelbrot-hue by Richard A Thomson
20. ↑phase by Linas Vepstas
21. ↑Complex numbers by David E Joyce
22. ↑Powers of complex numbers from Suitcase of Dreams
23. ↑Parabolic Julia Sets are Polynomial Time Computable by Mark Braverman
24. ↑SYMBOLIC DYNAMICS AND SELF-SIMILAR GROUPS by VOLODYMYR NEKRASHEVYCH
25. ↑ Query about general Julia set IFS for higher powers.
26. ↑wikipedia : repelling fixed point
27. ↑mathbbc/185430 mathoverflow question: clustering-of-periodic-points-for-a-polynomial-iteration-of-mathbbc
28. ↑Wolfram Alpha
29. ↑example
30. ↑wikipedia : repelling fixed point
31. ↑Fractint documentation - Inverse Julias
32. ↑Image and c source code : IIMM using hit limit
33. ↑Exploding the Dark Heart of Chaos by Chris King
34. ↑Drakopoulos V., Comparing rendering methods for Julia sets, Journal of WSCG 10 (2002), 155-161
35. ↑bugman123: Fractals
36. ↑dhushara : DarkHeart

• Book:Fractals