コンテンツにスキップ

英文维基 | 中文维基 | 日文维基 | 草榴社区

ファイル:Golden Mean Quadratic Siegel Disc Speed.png

ページのコンテンツが他言語でサポートされていません。

元のファイル (1,000 × 1,000 ピクセル、ファイルサイズ: 279キロバイト、MIME タイプ: image/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 )
日付
原典 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 .
作者 Adam majewski

ライセンス

この作品の著作権者である私は、この作品を以下のライセンスで提供します。
w:ja:クリエイティブ・コモンズ
表示 継承
このファイルはクリエイティブ・コモンズ 表示-継承 3.0 非移植ライセンスのもとに利用を許諾されています。
あなたは以下の条件に従う場合に限り、自由に
  • 共有 – 本作品を複製、頒布、展示、実演できます。
  • 再構成 – 二次的著作物を作成できます。
あなたの従うべき条件は以下の通りです。
  • 表示 – あなたは適切なクレジットを表示し、ライセンスへのリンクを提供し、変更があったらその旨を示さなければなりません。これらは合理的であればどのような方法で行っても構いませんが、許諾者があなたやあなたの利用行為を支持していると示唆するような方法は除きます。
  • 継承 – もしあなたがこの作品をリミックスしたり、改変したり、加工した場合には、あなたはあなたの貢献部分を元の作品とこれと同一または互換性があるライセンスの下に頒布しなければなりません。

Compare with

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

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 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

キャプション

このファイルの内容を1行で記述してください

このファイルに描写されている項目

題材

29 10 2011

ファイルの履歴

過去の版のファイルを表示するには、その版の日時をクリックしてください。

日付と時刻サムネイル寸法利用者コメント
現在の版2011年10月29日 (土) 12:222011年10月29日 (土) 12:22時点における版のサムネイル1,000 × 1,000 (279キロバイト)Soul windsurfer

以下のページがこのファイルを使用しています:

グローバルなファイル使用状況

以下に挙げる他のウィキがこの画像を使っています: