Fractals/Iterations in the complex plane/av velocity

Average velocity by Chris King "Discrete Velocity of non-attracting Basins and Petals: Compute, for the points that don't escape, the average discrete velocity on the orbit:" Combined Methods of Depicting Julia Sets by Chris KingExploding the Dark Heart of Chaos Chris King March 2009 – Jul 2016 Mathematics Department University of Auckland Algorithm right|thumb|distance between points 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 : 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 bandswikipedia : Colour banding 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 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. Images with src code ( see commons page) References