Non-analytic Fractional Iteration

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

Analytic Fractional Iteration

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

Real-analytic Iteration

Complex-analytic Iteration or Holomorphic Iteration

Forward iteration

The 10 first powers of a complex number inside the unit circle

Exponential spirals

Principle branch of arg

${\displaystyle z=re^{i\theta }\,}$

where :

• r is the absolute value or modulus or magnitude 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}}}.\,}$

${\displaystyle \theta =\arg(z)=\operatorname {atan2} (y,x)={\begin{cases}\arctan({\frac {y}{x}})&{\mbox{if }}x>0\\\arctan({\frac {y}{x}})+\pi &{\mbox{if }}x<0{\mbox{ and }}y\geq 0\\\arctan({\frac {y}{x}})-\pi &{\mbox{if }}x<0{\mbox{ and }}y<0\\{\frac {\pi }{2}}&{\mbox{if }}x=0{\mbox{ and }}y>0\\-{\frac {\pi }{2}}&{\mbox{if }}x=0{\mbox{ and }}y<0\\{\mbox{indeterminate }}&{\mbox{if }}x=0{\mbox{ and }}y=0.\end{cases}}}$

so

${\displaystyle f_{0}(z)=z^{2}=(re^{i\theta })^{2}=r^{2}e^{i2\theta }\,}$

and forward iteration :^[16]

${\displaystyle f_{0}^{n}(z)=r^{2^{n}}e^{i2^{n}\theta }\,}$

Forward iteration:

• squares radius and doubles angle ( phase, argument)^[17][18]
• gives forward orbit = list of points {z0, z1, z2, z3... , zn} which lays on exponential spirals.^[19][20]

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

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

${\displaystyle |z_{n}|>ER}$

Backward iteration

Backward iteration of complex quadratic polynomial with proper chose of the preimage

Every angle α ∈ R/Z measured in turns has :

• one image = 2α mod 1 under doubling map
• "two preimages under the doubling map: α/2 and (α + 1)/2.".^[22] 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 without extra informations.

Not that preimages show rotational symmetry ( 180 degrees)

For other functions see Fractalforum^[23]

Level curves of escape or attracting time

• 
  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

• 
  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^[24] of complex quadratic polynomial ${\displaystyle Z_{0}=\beta _{c}\,}$
• preimages of ${\displaystyle Z_{0}\,}$ by inverse 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[25]

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

See also :

• online tool : Wolfram Alpha^[26][27]

                                               beta
                    beta                                            -beta
   beta                         -beta                    x                     y

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

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

# octave m-file 
# IIM using Hit table 
# save as a iimm_hit.m in octave working dir
# run octave and iimm_hit
#
# stack-like operation
# http://www.gnu.org/software/octave/doc/interpreter/Miscellaneous-Techniques.html#Miscellaneous-Techniques
# octave can resize array 
# a = [];
# while (condition)
# a(end+1) = value; # "push" operation
# a(end) = []; # "pop" operation
#

# -------------- general octave settings ----------
clear all;
more off;
pkg load image; # imwrite
pkg load miscellaneous; # waitbar

# --------------- const and var -----------------------------
# define some global var AT EACH LEVEL where you want to use it ( Pascal CdeMills)
# ? for global variables one can't define and give initial value at the same time

 # integer ( screen ) coordinate
 iSide = 1000
 ixMax = iSide
 iyMax = iSide

# 
 global HitLimit;
 HitLimit=1 # proportional to number of details and time of drawing
 global Hits; # 2D array of pixels . Hit > 0 means that point was in orbit
 Hits = zeros(iyMax,ixMax); 
 
 
 # image as a 2D matrix of 24 bit colors coded from 0.0 to 1.0 
 global MyImage;
 MyImage = zeros(iyMax,ixMax,3); # matrix filled with 0.0 ( black image)= [rows, columns, color_depth]

# world ( float) coordinate - dynamical (Z) plane 
global  dSide;
global  ZxMin ;
global  ZxMax;
global ZyMin;
global ZyMax ;
global z;
global PixelHeight ;
global PixelWidth ;

# add values to global variables 
dSide=1.25
ZxMin = -dSide
ZxMax = dSide
ZyMin = -dSide
ZyMax = dSide
PixelHeight = (ZyMax - ZyMin)/(iyMax - 1)
PixelWidth = (ZxMax - ZxMin)/(ixMax - 1)

# pseudo stack = resizable array
global Stack;
Stack=[]; 
global StackIndex;
StackIndex=0;

c = complex(.27327401, 0.00756218); 
 
global Color24White;
Color24White=[1.0, 1.0, 1.0];

# ---------------- functions ------------------

function [iY, iX] = f(z)
  # define some global var AT EACH LEVEL where you want to use it ( Pascal CdeMills)
  global ZxMin;
  global ZyMax;
  global PixelWidth;
  global PixelHeight;

  # translate from world to screen coordinate
   iX=int32((real(z)-ZxMin)/PixelWidth);
   iY=int32((ZyMax - imag(z))/PixelHeight); # invert y axis

endfunction;

# plot point with integer coorfinate to array MyImage
function Plot( iY, iX , color )

  
  global MyImage;
  
   MyImage(iY, iX, 1) = color(1);
   MyImage(iY, iX, 2) = color(2);
   MyImage(iY, iX, 3) = color(3);
   
endfunction;

# push = put point z on the stack 
function push(z)

  global Stack;
  global StackIndex;
  
  StackIndex =size(Stack,2); # update global var 
  StackIndex = StackIndex +1; # update global var 
  Stack(StackIndex) =  z;  # "push" z on pseudo stack 
   
endfunction;

# pop = take point z from the stack
function z = pop()

 global Stack;
 global StackIndex;

 StackIndex =size(Stack,2); # update global var 
 z = Stack(StackIndex); # pop z from pseudo stack
 Stack(StackIndex) = []; # make pseudo stack shorter
 StackIndex = StackIndex -1 ; # update global var

endfunction;

function  CheckPoint(z)
# error is here 
 global Hits;
 global HitLimit;
 global Color24White;

 [iY, iX] = f(z);
 hit = Hits(iY, iX);
 if (hit < HitLimit)
   
   push(z); # (put z on the stack)  to continue iteration
   if (hit<1) Plot( iY, iX , Color24White ); endif; # plot 
   Hits(iY, iX) = hit + 1;   #  update Hits table
 endif;

endfunction; # CheckPoint

# -------------------- main ---------------------------------------

# Determine the unstable fixed point beta 
# http://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappings
beta = 0.5+sqrt(0.25-c)

z=-beta
CheckPoint(z);
while (StackIndex>0) # if stack is not empty 
 z = pop(); # take point z from the stack
 # compute its 2 preimages z and -z ( inverse function of complex quadratic polynomial)
 z= sqrt(z-c);
 # check points : draw, put on the stack to continue iterations 
 CheckPoint(z);
 CheckPoint(-z);
endwhile;

# ------------------- image --------------------------------------
 imread("a7.png");# first load any existing image to resolve bug : panic crash using imwrite: octave: magick/semaphore.c:525 first load any image
 image(MyImage); # Display Image 
 name =  strcat("iim",int2str(HitLimit), " .png");
 imwrite(MyImage,name);  # save image to file. this requires the ImageMagick "convert" utility.