Fractals/Iterations in the complex plane/boettcher

From Wikibooks, open books for an open world
Jump to: navigation, search

Intro[edit]

Superattracting fixed points[edit]

For complex quadratic polynomial there are many superattracting fixed point ( with multiplier = 0 ):

  • infinity ( It is always is superattracting fixed point for polynomials )
  • z_s = 0 \, is finite superattracting fixed point for map f_0\,
  • z_s = 0 \, and z_s = -1 \, are two finite superattracting fixed points for map f_{-1}\,

Description[edit]

Near[1] super attracting fixed point (for example infinity) the behaviour of discrete dynamical system :

z_{n+1} = f_c(z_n) = z_n^2 + c  \,

based on complex quadratic polynomial f_c(z) = z^2 + c\, is similar to

w_{n+1} = f_0(w_n) = w_n^2  \,

based on f_0(w) = w^2\,

It can be treated as one dynamical system viewed in two coordinate systems :

  • easy ( w )
  • hard to analyse( z )

[2]

In other words map f_c\, is conjugate [3] to map f_0\, near infinity. [4]

History[edit]

In 1904 LE Boettcher solved Schröder equation[5][6] in case of supperattracting fixed point[7]

Names[edit]

  • w\, is Boettcher coordinate
  • \Phi_c \, is Boettcher function
  • Boettcher Functional Equation [8][9]: \Phi_c(f_c(z)) = \Phi_c(z)^2 \,

where :

w = \Phi_c(z)\,

Exterior or complement of filled-in Julia set is :

A_{f_c}(\infty) = K(f_c)^C = F_\infty

It can be analysed using

  • escape time (simple but gives only radial values = escape time ) LSM/J,
  • distance estimation ( more advanced, continuus, but gives only radial values = distance ) DEM/J
  • Boettcher coordinate or complex potential ( the best , gives  :
    • radial values ( real potential ) LogPhi = CPM/J
    • angular values ( external angle ) ArgPhi

So both values can be used to color with 2D gradient.


First compute Boettcher coordinate w\, using this formula [10]

w = \Phi_c(z) = z* \prod_{n=1}^\infty \left( \frac{f_c^n(z)}{f_c^n(z)-c}\right )^{1/2^n}\,

It looks "simple", but square root of complex nuber gives two values so one have to choose one value.

Complex potential on dynamical plane[edit]

LogPhi - Douady-Hubbard potential - real potential - radial component of complex potential[edit]

CPM/J[edit]

Potential of filled Julia set
Diagram of potential computed with 2 methods : simple and full

Note that potential inside Kc is zero so :

Pseudocode version :

if (LastIteration==IterationMax)
 then potential=0    /* inside  Filled-in Julia set */
 else potential= GiveLogPhi(z0,c,ER,nMax); /* outside */

It also removes potential error for log(0).

Full version[edit]

Math (full) notation : [11]

LogPhi_c(z) =  ln |z| + \sum_{n=1}^\infty \frac{1}{2^n} ln | 1 + \frac{c}{(f_c^{n-1}(z))^2}|


Maxima (full) function :

GiveLogPhi(z0,c,ER,nMax):=
block(
 [z:z0,
 logphi:log(cabs(z)),
 fac:1/2,
 n:0],
 while n<nMax and abs(z)<ER do
 (z:z*z+c,
 logphi:logphi+fac*log(cabs(1+c/(z*z))),
 n:n+1
 ),
 return(float(logphi))
)$

Simplified version[edit]

The escape rate function of a polynomial f is defined by :

G_f(z) = \lim_{n\rightarrow \infty}\frac{1}{2^n} log^+ |f^n(z)| \,

where :

log^+ = max(log,0)

"The function Gp is continous on C and harmonic on the complement of the Julia set. It vanishes identically on K(f) and as it has a logarithmic pole at infinity, it is a it is the Green's function for C/ K(f)." ( Laura G. DeMarco) [12]



Math simplified formula :

SLogPhi_c(z) = \frac{log(f^n(z))}{ 2^n} \,


Maxima function :

GiveSLogPhi(z0,c,e_r,i_max):=
 block(
  [z:z0,
   logphi,
   fac:1/2,
   i:0
  ],
  while i<i_max and cabs(z)<e_r do
  (z:z*z+c,
    fac:fac/2,
    i:i+1
  ),
 logphi:fac*log(cabs(z)),
 return(float(logphi))
)$

If you don't check if orbit is not bounded ( escapes, bailout test) then use this Maxima function :

GiveSLogPhi(z0,c,e_r,i_max):=
block(
 [z:z0, logphi, fac:1/2, i:0],
  while i<i_max and cabs(z)<e_r do
   (z:z*z+c,
    fac:fac/2,
    i:i+1 ),
   if i=i_max
    then logphi:0
    else logphi:fac*log(cabs(z)),
   float(logphi) 
)$


C version :


double jlogphi(double zx0, double zy0, double cx, double cy)
/* this function is based on function by W Jung http://mndynamics.com */
{ 
 int j; 
 double 
 zx=zx0,
 zy=zy0,
 s = 0.5, 
 zx2=zx*zx,
 zy2=zy*zy,
 t;
 for (j = 1; j < 400; j++)
 { s *= 0.5; 
  zy = 2 * zx * zy + cy;
  zx = zx2 - zy2 + cx; 
  zx2 = zx*zx; 
  zy2 = zy*zy;
  t = fabs(zx2 + zy2); // abs(z)
  if ( t > 1e24) break; 
 } 
return s*log2(t);  // log(zn)* 2^(-n)
}//jlogphi


Euler version by R. Grothmann ( with small change : from z^2-c to z^2+c) :[13]

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

x=-2:0.05:2; y=x'; z=x+I*y;
{w,n}=iter(z,c);
wr=max(0,log(abs(w)))/2^n;

Level Sets of potential = pLSM/J[edit]

Here is Delphi function which gives level of potential :

Function GiveLevelOfPotential(potential:extended):integer;
 var r:extended;
 begin
    r:= log2(abs(potential));
    result:=ceil(r);
 end;

Level Curves of potential = equipotential lines = pLCM/J[edit]

ArgPhi - External angle - angular component of complex potential[edit]

One can start with binary decomposition of basin of attraction of infinity.

The second step can be using \Psi_c\,

polar coordinate system and Psi_c for c=-2

period detection[edit]

How to find period of external angle measured in turns under doubling map :

Here is Common Lisp code :

(defun give-period (ratio-angle)
  "gives period of angle in turns (ratio) under doubling map"
  (let* ((n (numerator ratio-angle))
	 (d (denominator ratio-angle))
	 (temp n)) ; temporary numerator
 
    (loop for p from 1 to 100 do 
	  (setq temp  (mod (* temp 2) d)) ; (2 x n) modulo d = doubling)
	  when ( or (= temp n) (= temp 0)) return p )))


Maxima CAS code :

doubling_map(n,d):=mod(2*n,d);

/* catch-throw version by Stavros Macrakis, works */
GivePeriodOfAngle(n0,d):=
catch(
      block([ni:n0],
          for i thru 200 do if (ni:doubling_map(ni,d))=n0 then throw(i),
          0 ) )$


/* go-loop version, works */
GiveP(n0,d):=block(
[ni:n0,i:0],
block(
  loop,
  ni:doubling_map(ni,d),
  i:i+1,
  if i<100 and not (n0=ni) then go(loop)
),
if (n0=ni) 
        then i 
        else 0
);

/* Barton Willis while version without for loop , works */
GivePeriod(n0,d):=block([ni : n0,k : 1],
  while (ni : doubling_map(ni,d)) # n0 and k < 100 do (
    k : k + 1),
  if k = 100 then 0 else k)$


Computing external angle

External angle (argument) is argument of Boettcher coordinate w\,

arg_c(z) = arg(w) = arg(\Phi_c(z)) \,

Because Boettcher coordinate is a product of complex numbers so argument of product is :


arg_c(z) = arg(z) + \sum_{n=1}^\infty \left( \frac{1}{2^n}*arg \left( \frac{f_c^n(z)}{f_c^n(z)-c}     \right ) \right )  \,


Constructing the spine of filled Julia set[edit]

Algorithm for constructiong the spine is described by A. Douady[14]

  • join - \beta\, and \beta\,,
  • (to do )

Drawing dynamic external ray[edit]

Field lines in in the Fatou domain[edit]

Explanation by Gert Buschmann

backwards iteration[edit]

Backward iteration of complex quadratic polynomial with proper choice of the preimage
Periodic external rays landing on alfa fixed points for periods 2-40. Made with backwards iteration

This method has been used by several people and proved by Thierry Bousch. [15]

Code in c++ by Wolf Jung can be found in procedure QmnPlane::backray() in file qmnplane.cpp ( see source code of program mandel version 5.3 ). [16]


  • Ray for periodic angle ( simplest case )

It will be explained by an example :

First choose external angle t \, (in turns). External angle for periodic ray is a rational number.

t = \frac{1}{3}\,

Compute period of external angle under doubling map.

Because "1/3 doubled gives 2/3 and 2/3 doubled gives 4/3, which is congruent to 1/3" [17]

1 \equiv 2 \pmod 3.\,

or

\frac{2}{3} \equiv \frac{1}{3} \pmod 1.\, [18]

so external angle t = \frac{1}{3}\, has period 2 under doubling map.

Start with 2 points near infinity (in conjugate plane):

on ray 1/3 is a point w = R*e^{2\pi i/3}\,

on ray 2/3 is a point w = R*e^{4 \pi i /3}\,.

Near infinity z = w \, so one can swith to dynamical plane ( Boettcher conjugation )

Backward iteration (with proper chose from two possibilities)[19] of point on ray 1/3 goes to ray 2/3, back to 1/3 and so on.

In C it is :

/* choose one of 2 roots: zNm1 or -zNm1  where zN = sqrt(zN - c )  */
if (creal(zNm1)*creal(zN) + cimag(zNm1)*cimag(zN) <= 0) zNm1=-zNm1;

or in Maxima CAS :

if (z1m1.z01>0) then z11:z1m1 else z11:-z1m1;


One has to divide set of points into 2 subsets ( 2 rays). Draw one of these 2 sets and join the points. It will be an approximation of ray.

  • Ray for preperiodic angle ( to do )

Drawing dynamic external ray using inverse Boettcher map by Curtis McMullen[edit]

Julia set with external rays using McMullen method

This method is based on C program by Curtis McMullen[20] and its Pascal version by Matjaz Erat[21]

It consist of 3 big steps :

  • compute some w-points of external ray of circle for angle \theta  \, and various radii (rasterisation)
w_0 = R * e^{2\pi i t} \, where 1 < R < \infty \,
  • map w-points to z-point using inverse Boettcher map \Psi_c = \Phi_{c}^{-1} \,
z_0 = \Psi_c (w_0) = f_c^{-n}(w_0^{2^n}) \,
  • draw z-points ( and connect them using lines )


First and last steps are easy, but second is not so needs more explanation.

Rasterisation[edit]

For given external ray in f_0\, plane each point of ray has :

  • constant value t\, ( external angle in turns )
  • variable radius R\,

so w\, points of ray are parametrised by radius R\, and can be computed using exponential form of complex numbers :

w = F_t(R) = R * e^{2*\Pi*i*t}\,

One can go along ray using linear scale :

t:1/3; /* example value */
R_Max:4;
R_Min:1.1;
for R:R_Max step -0.5 thru R_Min do w:R*exp(2*%pi*%i*t);
/* Maxima allows non-integer values in for statement */ 

It gives some w points with equal distance between them.


Another method is to use nonlinera scale.

To do it we introduce floating point exponent r\, such that :

R = 2^r\,

and

w = G_t(r) = 2^r * e^{2*\Pi*i*t}\,

To compute some w points of external ray in f_0\, plane for angle t\, use such Maxima code :

t:1/3; /* external angle in turns */
/* range for computing  R ; as r tends to 0 R tends to 1 */
rMax:2; /* so Rmax=2^2=4 /
rMin:0.1; /* rMin > 0   */
caution:0.93; /* positive number < 1 ; r:r*caution  gives smaller r */
r:rMax;
unless r<rMin do
( 
 r:r*caution, /* new smaller r */ 
 R:2^r, /* new smaller R */
 w:R*exp(2*%pi*%i*t) /* new point w in f0 plane */
);

In this method distance between points is not equal but inversely proportional to distance to boundary of filled Julia set.

It is good because here ray has greater curvature so curve will be more smooth.

Mapping[edit]

Mapping points from f_0\,-plane to f_c\,-plane consist of 4 minor steps :

  • forward iteration in f_0\, plane
w_n = w_0^{2^n} = R^{2^n}*e^{2\pi i 2^n t}\,

until w_n\, is near infinity

abs(w_n) > R_{limit} \,
  • switching plane ( from f_0\, to f_c\,)
z_n =w_n \,

( because here, near infinity : w = z \,)

  • backward iteration in f_c\, plane the same ( n \,) number of iterations
  • last point  z_0 \, is on our external ray


1,2 and 4 minor steps are easy. Third is not.

z_{n-1} = f_c^{-1} (z_n) = \sqrt{z_n - c}\,

Backward iteration uses square root of complex number. It is 2-valued functions so backward iteration gives binary tree.

One can't choose good path in such tree without extre informations. To solve it we will use 2 things :

  • equicontinuity of basin of attraction of infinity
  • conjugacy between f_0\, and f_c\, planes
Equicontinuity of basin of attraction of infinity[edit]

Basin of attraction of infinity ( complement of filled-in Julia set) contains all points which tends to infinity under forward iteration.

A_{f_c}(\infty) \  \overset{\underset{\mathrm{def}}{}}{=} \  \{ z \in  \mathbb{C}  : f^{(k)} _c (z)  \to  \infty\  as\  k \to \infty \}.

Infinity is superattracting fixed point and orbits of all points have similar behaviour. In other words orbits of 2 points are assumed to stay close if they are close at the beginning.

It is equicontinuity ( compare with normality).

In f_c\, plane one can use forward orbit of previous point of ray for computing backward orbit of next point.

Detailed version of algorithm[edit]
  • compute first point of ray (start near infinity ang go toward Julia set )
w = 2^r * e^{2\Pi i t}\, where r = r_{max}\,

here one can easily switch planes :

z = w\,

It is our first z-point of ray.

  • compute next z-point of ray
    • compute next w-point of ray for r = r * caution \,
    • compute forward iteration of 2 points : previous z-point and actual w-point. Save z-orbit and last w-point
    • switch planes and use last w-point as a starting point :z_n = w_{last}\,
    • backward iteration of new z_n\, toward new z_0\, using forward orbit of previous z point
    • z_0\, is our next z point of our ray
  • and so on ( next points ) until  r < r_{Min}\,

Maxima CAS src code

/* gives a list of z-points of external ray for angle t in turns and coefficient c */
GiveRay(t,c):=
block(
 [r],
 /* range for drawing  R=2^r ; as r tends to 0 R tends to 1 */
 rMin:1E-20, /* 1E-4;  rMin > 0  ; if rMin=0 then program has infinity loop !!!!! */
 rMax:2,  
 caution:0.9330329915368074, /* r:r*caution ; it gives smaller r */
 /* upper limit for iteration */
 R_max:300,
 /* */
 zz:[], /* array for z points of ray in fc plane */
 /*  some w-points of external ray in f0 plane  */
 r:rMax,
 while 2^r<R_max do r:2*r, /* find point w on ray near infinity (R>=R_max) in f0 plane */
 R:2^r,
 w:rectform(ev(R*exp(2*%pi*%i*t))),
 z:w, /* near infinity z=w */
 zz:cons(z,zz), 
 unless r<rMin do
 ( /* new smaller R */
  r:r*caution,  
  R:2^r,
  /* */
  w:rectform(ev(R*exp(2*%pi*%i*t))),
  /* */
  last_z:z,
  z:Psi_n(r,t,last_z,R_max), /* z=Psi_n(w) */
  zz:cons(z,zz)
 ),
 return(zz)
)$

Lamination of dynamical plane[edit]

Lamination of rabbit Julia set

Here is long description

See also[edit]

References[edit]

  1. Neighbourhood in wikipedia
  2. The work of George Szekeres on functional equations by Keith Briggs
  3. Topological conjugacy in wikipedia
  4. How to draw external rays by Wolf Jung
  5. Schröder equation in wikipedia
  6. Lucjan Emil Böttcher and his mathematical legacy by Stanislaw Domoradzki, Malgorzata Stawiska
  7. L. E. Boettcher, The principal laws of convergence of iterates and their aplication to analysis (Russian), Izv. Kazan. fiz.-Mat. Obshch. 14) (1904), 155-234.
  8. Böttcher equation at Hyperoperations Wiki
  9. wikipedia : Böttcher's equation
  10. How to draw external rays by Wolf Jung
  11. The Beauty of Fractals, page 65
  12. Holomorphic families of rational maps: dynamics, geometry, and potential theory. A thesis presented by Laura G. DeMarco
  13. Euler examples by R. Grothmann
  14. A. Douady, “Algorithms for computing angles in the Mandelbrot set,” in Chaotic Dynamics and Fractals, M. Barnsley and S. G. Demko, Eds., vol. 2 of Notes and Reports in Mathematics in Science and Engineering, pp. 155–168, Academic Press, Atlanta, Ga, USA, 1986.
  15. Thierry Bousch : De combien tournent les rayons externes? Manuscrit non publié, 1995
  16. Program Mandel by Wolf Jung
  17. Explanation by Wolf Jung
  18. Modular arithmetic in wikipedia
  19. Square root of complex number gives 2 values so one has to choose only one. For details see Wolf Jung page
  20. c program by Curtis McMullen (quad.c in Julia.tar.gz)
  21. Quadratische Polynome by Matjaz Erat