Summary

Licensing

I, the copyright holder of this work, hereby publish it under the following license:

This file is licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license.

You are free:

• to share – to copy, distribute and transmit the work
• to remix – to adapt the work

Under the following conditions:

• attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
• share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.

https://creativecommons.org/licenses/by-sa/3.0CC BY-SA 3.0 Creative Commons Attribution-Share Alike 3.0 truetrue

Compare with

• 
  Gray and detailed version
• 
  Boundary made with MIIM
• 
  Animated version
• 
  Orbits inside Siegel Disc

Theory

Program draws dynamic plane for discrete dynamical system based on complex quadratic polynomial^[3]:

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

Rotation number ( internal angle) t is an irrational number = the Golden Mean :^[4]

${\displaystyle t={\frac {{\sqrt {5}}-1}{2}}\approx 0.618033988749895}$

It is used to compute c on the boundary of main cardioid  :^[5]

${\displaystyle c=c(t)={\frac {e^{\Pi *t*i}}{2}}-{\frac {e^{2*\Pi *t*i}}{4}}\approx -0.390540870218399-0.586787907346969*i}$

Inside of filled Julia set is a Siegel Disc^[6] around fixed point alpha^[7]:

${\displaystyle z=\alpha _{c}={\frac {1-{\sqrt {1-4c}}}{2}}=-0.36868-0.33775i}$   

with multiplier

${\displaystyle \lambda (\alpha _{c})=1-{\sqrt {1-4c}}\,}$

such that

${\displaystyle abs(\lambda (\alpha _{c}))=1.0000\,}$

Algorithm

On dynamical plane one can see :

• Exterior of filled Julia set (blue) colored by level set method,
• Interior of Julia set showing irrational flow (green) coloured by the sine of the velocity

For the points that don’t escape compute the average discrete velocity of orbit  :

${\displaystyle \upsilon ={\frac {S_{n}}{n}}}$

where :

${\displaystyle S_{n}=\sum _{n=1}^{n_{max}}d_{n}=\sum _{n=1}^{n_{max}}|Z_{n+1}-Z_{n}|}$

In Octave it looks :

# octave code
 d=0;
 iter = 0;

 while (iter < maxiter) && (abs(z)<ER)
   h=z; # previous point = z_(n)
   z=z*z+c; # next point = z_(n+1)
   iter = iter+1;
   d=d+abs(z-h); # sum of distances along orbit
 end

 if iter < maxiter  # exterior 
    measure = iter;
    myflag=3; # escaping to infinity 
             
    else # iter==maxiter ( inside filled julia set )
      measure=20*d/iter; # average distance (d/iter) = 0.5

In Chris King Maple code this discrete velocity is measured only by sum of distances between points

${\displaystyle S_{n}=\sum _{n=1}^{n_{max}}d_{n}=\sum _{n=1}^{n_{max}}|Z_{n+1}-Z_{n}|}$

Because :

• all forward orbit from interior of Julia set fall into SIegel disc
• inside Siegel disc points turn around its center ( indifferent periodic point )

so distance is a good measure into which Siegel orbit point fall

Using periodic function ( sin, cos) creates bands^[8] showing dynamics inside Julia set ( siegel disc and its preimages ).

Matlab src code

% code by Chris King
% http://www.dhushara.com/DarkHeart/Viewers/source/siegel.m
function siegel();
nx = 480;
ny = 480;
ColorMset = zeros(nx,ny,3);
magc=0.65;
xmin = -1/magc;
xmax = 1/magc;
ymin = -1/magc;
ymax = 1/magc;
maxiter = 1200;
wb = waitbar(0,'Please wait...');
for iy = 1:ny
  cy = ymin + iy*(ymax - ymin)/(ny - 1);
  for ix= 1:nx
    cx = xmin + ix*(xmax - xmin)/(nx - 1);
    [k myfl] = Mlevel(cy,cx,maxiter);
    if myfl==2  
        ColorMset(ix,iy,2) = abs(sin(5*k/10+pi/4));
    else
        if myfl==1
            ColorMset(ix,iy,1) = abs(sin(2*k/10));
        else
            %ColorMset(ix,iy,2) = abs(sin(2*k/10+pi/4));
            ColorMset(ix,iy,3) = abs(cos(2*k/10));
        end
    end
  end
  waitbar(iy/ny,wb)
end
close(wb);
image(ColorMset);
imwrite(ColorMset,'siegel.jpg','jpg','Quality',100);
 
function [potential myfl] = Mlevel(cx,cy,maxiter)
z = complex(cx,cy);
th=pi*(-1+sqrt(5));
d=exp(complex(0,th));
d=d/2-d*d/4;
%e=(1-sqrt(1-4*d))/2;
%e=0;
%a=complex(0,sqrt(3));
%a=sqrt(3);
a=4;
ang=0;
iter = 0;
while (iter < maxiter)&&(abs(z) > 0.001)&&(abs(z)<20)
   h=z;
   %z=d*z*z*(z-a)/(1-a*z);
   z=z*z+d;
     hh=abs(z-h)*(z-h);
     if iter>maxiter/2
       ang=ang+hh;
    end
   iter = iter+1;
end
if iter < maxiter
    potential = iter;
    if abs(z)>=20
        myfl=0;
    else
        myfl=1;
    end
else
    %potential = -(ang-floor(ang));
    potential=abs(ang);
    myfl=2;
end

Changes / questions

1. orientation of the plane : In Matlab code is definition :

function [potential myfl] = Mlevel(cx,cy,maxiter)

but use is different :

[k myfl] = Mlevel(cy,cx,maxiter);

I have changed :

[k myfl] = Mlevel(cx,cy,maxiter); # order of arguments
ColorMset = zeros(ny,nx,3); # order of arguments
cy = ymax - iy*(ymax - ymin)/(ny - 1); # reverse y axis
ColorMset(iy,ix,1) = abs(sin(2*k/10)) # order of arguments

Now the orientation is the same as in this image^[9]

I check it with :

if(cy>0 && cx>0)  ColorMset(iy,ix,2)=1.0-MyImage(iy,ix,2);
 

2. Output file type : Png file is better then jpg in case of raster graphic

3. Speed

Is it posible to vectorize computations like in this r code^[10]?

4. Names of variables

Octave src code

# http://www.dhushara.com/DarkHeart/DarkHeart.htm
# it is Octave m-file
# converted from matlab m-file by Chris King
# http://www.dhushara.com/DarkHeart/Viewers/source/siegel.m
#

# ------------- load packages ------------------------
pkg load image;
pkg load miscellaneous; # waitbar

# --------- definitions ------------------------------
 
function [potential myfl] = Mlevel(zx,zy,c,maxiter)
 ER=2.0; # escape radius = bailout value 
 z = complex(zx,zy);
 ang=0;
 iter = 0;

 while (iter < maxiter) && (abs(z) > 0.001) && (abs(z)<ER)
   h=z; # previous point = z_(n)
   z=z*z+c; # next point = z_(n+1)
   

     # for the points that don''t escape compute 
     #  the average discrete velocity  on the orbit = abs( z_(n+1) - z_n ) 
     if iter>maxiter/2 # ???
       zh=z-h; 
       hh=abs(zh)*zh;
       ang=ang+hh;
    endif;

   iter = iter+1;
 end

 if iter < maxiter 
    potential = iter;
    if abs(z)>=ER  myfl=3; # escaping to infinity 
             else  myfl=1; # ??? falling into Siegel disc
    end
 else # iter==maxite ( inside filled julia set )
    potential=abs(ang);
    myfl=2;
 end
endfunction; # Mlevel

# ------------- const ------------------------------

# integer ( screen ) coordinate 
iSide=1000
nx = iSide;
ny = iSide;

# image as a 2D matrix of 24 bit colors coded from 0.0 to 1.0 
MyImage = zeros(ny,nx,3); # matrix filled with 0

# world ( float) coordinate - dynamical (Z) plane 
magc=0.65;
dSide=1/magc
Zxmin = -dSide;
Zxmax = dSide;
Zymin = -dSide;
Zymax = dSide;

stepy = (Zymax - Zymin)/(ny - 1); # pixel height
stepx = (Zxmax - Zxmin)/(nx - 1); # pixel width

maxiter = 2000

# fc(z) = z*z + c 
# rotation number or internal angle
t = (-1+sqrt(5))/2
th=2*pi*t;  # from turns to radians
d=exp(complex(0,th)); # 
c =d/2-d*d/4 # point on boundary of main cardioid


pi4=pi/4;

# --------------- main : 2 loops ---------------------

waitbar(0,'Please wait...'); # info 

# scan all pixels of image and comput color 
for iy = 1:ny
  Zy = Zymax - iy*stepy; # invert y axis
  for ix= 1:nx
    Zx = Zxmin + ix*stepx; # map from screen to world coordinate

    [k myfl] = Mlevel(Zx,Zy,c, maxiter);
    # color 
    switch (myfl)
     case 1 MyImage(iy,ix,1) = abs(sin(2*k/10)); # ?? Julia set in red
     case 2 MyImage(iy,ix,2) = abs(sin(5*k/10 + pi4)); # irrational flow (green) by the sine of the velocity.
     case 3 MyImage(iy,ix,3) = abs(cos(2*k/10)); # Exterior (blue) by level sets of escape time
    endswitch;
   
   # check plane orientation
   # first quadrant should be in upper right position
   # if(Zy>0 && Zx>0)  
   # MyImage(iy,ix,2)=1.0-MyImage(iy,ix,2);
   # endif;
  endfor; # for ix
  
  waitbar(iy/ny);
endfor; # for iy

# image 

image(MyImage); # display image
imwrite(MyImage,'si-test.png');  # save image to file
# this requires the ImageMagick "convert" utility.

References

1. ↑ Image of Julia sets by Chris King
2. ↑ Matlab m-file siegel.m by Chris King
3. ↑ wikipedia : Complex quadratic polynomial
4. ↑ Golden ratio at wikipedia
5. ↑ wikibooks : boundary_of_main_cardioid
6. ↑ siegel disc at wikibooks
7. ↑ Periodic points of complex quadratic mappings at wikipedia
8. ↑ wikipedia : Colour banding
9. ↑ Siegel disc image
10. ↑ Mandelbrot animation in R