Original file (1,000 × 1,000 pixels, file size: 279 KB, MIME type: image/png)
This is a file from the Wikimedia Commons. Information from its description page there is shown below. Commons is a freely licensed media file repository. You can help. |
Contents
Summary
DescriptionGolden Mean Quadratic Siegel Disc Speed.png |
English: Golden Mean Quadratic Siegel Disc with interior colored propotional to the average discrete velocity on the orbit = abs( z_(n+1) - z_n ) |
Date | |
Source | It is a copy of image[1] by Chris King made with use of his Matlab code[2]. ( note that Chris have made also Julia set (white) by modified inverse iteration ). I have only converted code to Octave and made some small changes. The core of algorithm remains unchanged. Thx to Chris for the great code and releasing it under free licence . |
Author | Adam majewski |
Licensing
- 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.
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]:
Rotation number ( internal angle) t is an irrational number = the Golden Mean :[4]
It is used to compute c on the boundary of main cardioid :[5]
Inside of filled Julia set is a Siegel Disc[6] around fixed point alpha[7]:
with multiplier
such that
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 :
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[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
- ↑ Image of Julia sets by Chris King
- ↑ Matlab m-file siegel.m by Chris King
- ↑ wikipedia : Complex quadratic polynomial
- ↑ Golden ratio at wikipedia
- ↑ wikibooks : boundary_of_main_cardioid
- ↑ siegel disc at wikibooks
- ↑ Periodic points of complex quadratic mappings at wikipedia
- ↑ wikipedia : Colour banding
- ↑ Siegel disc image
- ↑ Mandelbrot animation in R
Items portrayed in this file
depicts
some value
29 October 2011
File history
Click on a date/time to view the file as it appeared at that time.
Date/Time | Thumbnail | Dimensions | User | Comment | |
---|---|---|---|---|---|
current | 12:22, 29 October 2011 | 1,000 × 1,000 (279 KB) | Soul windsurfer |
File usage
The following page uses this file:
Global file usage
The following other wikis use this file:
- Usage on en.wikibooks.org
- Usage on ja.wikipedia.org