Fractals/Iterations in the complex plane/av velocity
< Fractals
Average velocity by Chris King
edit"Discrete Velocity of non-attracting Basins and Petals: Compute, for the points that don't escape, the average discrete velocity on the orbit:" [1][2]
Algorithm
editOn 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 :
where :
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
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[3] showing dynamics inside Julia set ( siegel disc and its preimages ).
Matlab src code
edit% 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
Octave src code
edit# 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.
Images with src code ( see commons page)
edit-
gray Siegel Disk
-
color Siegel DIsk