From maxima

• File
• File history
• File links

Julia set and external rays landing on parabolic period 2 cycle

Contents 1 Algorithm 1.1 How to compute coefficient ? 2 How to run a program 3 Maxima source code

Algorithm

Program draws (to png file) :

• dynamic rays using
• Julia set ( backward orbit of repeller = IIM/J )

How to compute coefficient ?

We want point c for which period 2 cycle z₀,z₁ is parabolic ( indifferent ). It is when absolute value of multiplier at periodic point is equal to one :

First compute point on unit circle with internal angle ( any other rational angle will also be good) :

Note that radius=1 gives abs(multiplier)=1

Map point to point on boundary of period 2 hyperbolic component of Mandelbrot set ( parameter plane ) using inverse multiplier map :

so

Our c is a root point of period 6 component connected with period 2 component.

How to run a program

From console under windows run (example commands):

cd C:\Program Files\Maxima-5.16.3\bin
maxima
batch("D:/t.mac")$

Maxima source code

This is batch file for Maxima CAS

/* --------------------------definitions of functions ------------------------------*/
f(z,c):=z*z+c;
finverseplus(z,c):=sqrt(z-c);
finverseminus(z,c):=-sqrt(z-c);
/*
Square root of complex number : csqrt(x + y * i) = sqrt((r + x) / 2) + i * y / sqrt(2 * (r + x))
gives principal value of square root 
*/
csqrt(z):=
 block(
 [t,re,im],
 t:abs(z)+realpart(z),
 if t>0 
  then (re:sqrt(t/2), im:imagpart(z)/sqrt(2*t))
  else  (im:abs(z), re:0),
 return(float(re+im*%i))
 )$
/*   */
Psi_n(r,t,z_last, Max_R):=
block(
 [iMax:200,
 iMax2:0],
 /* -----  forward iteration of 2 points : z_last and w --------------*/
 array(forward,iMax-1), /* forward orbit of z_last for comparison */
 forward[0]:z_last,
 i:0,
 while cabs(forward[i])<Max_R  and  i< ( iMax-2) do
 (	
  /* forward iteration of z in fc plane & save it to forward array */
  forward[i+1]:forward[i]*forward[i] + c, /* z*z+c */
  /* forward iteration of w in f0 plane :  w(n+1):=wn^2 */
  r:r*2, /* square radius = R^2=2^(2*r) because R=2^r */
  t:mod(2*t,1),
  /* */
 iMax2:iMax2+1,
 i:i+1
 ),
/* compute last w point ; it is equal to z-point */
R:2^r,
/* w:R*exp(2*%pi*%i*t),	z:w, */
array(backward,iMax-1),
backward[iMax2]:rectform(ev(R*exp(2*%pi*%i*t))), /* use last w as a starting point for backward iteration to new z */
/* -----  backward iteration point  z=w in fc plane --------------*/
for i:iMax2 step -1 thru 1 do
 (
  temp:csqrt(backward[i]-c), /* sqrt(z-c) */
  scalar_product:realpart(temp)*realpart(forward[i-1])+imagpart(temp)*imagpart(forward[i-1]),
  if (0>scalar_product) then temp:-temp, /* choose preimage */
   backward[i-1]:temp
  ),
return(backward[0])
)$
/*  */
GiveRay(t,c):=
block(
[r],
/* range for drawing  R=2^r ; as r tends to 0 R tends to 1 */
rMin:1E-20, /* 1E-4;  rMin > 0  ; if rMin=0 then program has infinity loop !!!!! */
rMax:2,  
caution:0.9330329915368074, /* r:r*caution ; it gives smaller r */
/* upper limit for iteration */
R_max:300,
/* */
zz:[], /* array for z points of ray in fc plane */
/*  some w-points of external ray in f0 plane  */
r:rMax,
while 2^r<R_max do r:2*r, /* find point w on ray near infinity (R>=R_max) in f0 plane */
R:2^r,
w:rectform(ev(R*exp(2*%pi*%i*t))),
z:w, /* near infinity z=w */
zz:cons(z,zz), 
unless r<rMin do
( /* new smaller R */
 r:r*caution,  
 R:2^r,
 /* */
 w:rectform(ev(R*exp(2*%pi*%i*t))),
 /* */
 last_z:z,
 z:Psi_n(r,t,last_z,R_max), /* z=Psi_n(w) */
 zz:cons(z,zz)
 ),
 return(zz)
)$
/* Gives points of backward orbit of z=repellor       */
GiveBackwardOrbit(c,repellor,zxMin,zxMax,zyMin,zyMax,iXmax,iYmax):=
block(
hit_limit:10, /* proportional to number of details and time of drawing */
PixelWidth:(zxMax-zxMin)/iXmax,
PixelHeight:(zyMax-zyMin)/iYmax,
/* 2D array of hits pixels . Hit > 0 means that point was in orbit */
array(Hits,fixnum,iXmax,iYmax), /* no hits for beginning */
/* choose repeller z=repellor as a starting point */
stack:[repellor], /*save repellor in stack */
/* save first point to list of pixels  */ 
x_y:[repellor], 
/* reversed iteration of repellor */
loop,
/* pop = take one point from the stack */
z:last(stack),
stack:delete(z,stack),
/*inverse iteration - first preimage (root) */
z:finverseplus(z,c),
/* translate from world to screen coordinate */
iX:fix((realpart(z)-zxMin)/PixelWidth),
iY:fix((imagpart(z)-zyMin)/PixelHeight),
hit:Hits[iX,iY],
if hit<hit_limit   
 then 
 (
 Hits[iX,iY]:hit+1,
 stack:endcons(z,stack), /* push = add z at the end of list stack */
 if hit=0 then x_y:endcons( z,x_y)
 ),
/*inverse iteration - second preimage (root) */
z:-z,
/* translate from world to screen coordinate, coversion to integer */
iX:fix((realpart(z)-zxMin)/PixelWidth),
iY:fix((imagpart(z)-zyMin)/PixelHeight),
hit:Hits[iX,iY],
if hit<hit_limit   
 then 
  (
   Hits[iX,iY]:hit+1,
   stack:endcons(z,stack), /* push = add z at the end of list stack to continue iteration */
   if hit=0 then x_y:endcons( z,x_y)
  ),
if is(not emptyp(stack)) then go(loop), 
return(x_y) /* list of pixels in the form [z1,z2] */
)$
/* ----------------------- main ----------------------------------------------------*/
start:elapsed_run_time ();
/* c:-0.12256+0.74486*%i;  value by Milnor*/
/* c:0.74486176661974*%i-0.12256116687665;  center of period 3 component */
c:float(rectform( exp(2*%pi*%i/3)/4 -1)); /* root point of period 6 component connected with period 2 component 0.21650635094611*%i-1.125*/
/* resolution is proportional to number of details and time of drawing */
iX_max:1000;
iY_max:1000;
/* define z-plane ( dynamical ) */
ZxMin:-2.0;
ZxMax:2.0;
ZyMin:-2.0;
ZyMax:2.0;
/* compute ray points & save to zz lists */
zz1:GiveRay(11/63,c); 
zz2:GiveRay(22/63,c);
zz3:GiveRay(25/63,c);  
zz4:GiveRay(37/63,c); 
zz5:GiveRay(44/63,c);
zz6:GiveRay(50/63,c); 
/* period 2 parabolic cycle, computed in other batch file */
z0:0.1703125096583*%i-1.135614939209353;
z1:0.13561493920935-0.1703125096583*%i;
/* save it to list zp */
zp:[];
zp:cons(z0,zp);
zp:cons(z1,zp);
beta:rectform((1+csqrt(1-4*c))/2); /* compute repelling fixed point beta */
alfa:rectform((1-csqrt(1-4*c))/2); /* other fixed point */
/* compute backward orbit of repelling fixed point 
xy: GiveBackwardOrbit(c,beta,ZxMin,ZxMax,ZyMin,ZyMax,iX_max,iY_max); /**/
*/
/* time of computations */
time:fix(elapsed_run_time ()-start);
/* draw it using draw package by */
load(draw); 
draw2d(
 terminal  = 'png,
 file_name = "JuliaRay_2_5",
 user_preamble="set size square;set key bottom right",
 title= concat("Dynamical plane for fc(z)=z*z+",string(c),"; Julia set and external
 rays landing on  period 2 parabolic orbit"),
 pic_width  = iX_max,
 pic_height = iY_max,
 yrange = [ZyMin,ZyMax],
 xrange = [ZxMin,ZyMax],
 xlabel     = "Z.re ",
 ylabel     = "Z.im",
 point_type = filled_circle,
 point_size    =  0.2,
 points_joined =false,
 color         = black,
 key = "backward orbit of z=beta",
 points(map(realpart,xy),map(imagpart,xy)),
 points_joined =true,
 color         = red,
 key = concat("external ray for angle ",string(11/63)),
 points(map(realpart,zz1),map(imagpart,zz1)),
 key = concat("external ray for angle ",string(22/63)),
 points(map(realpart,zz2),map(imagpart,zz2)),
 key = concat("external ray for angle ",string(25/63)),
 points(map(realpart,zz3),map(imagpart,zz3)),
 key = concat("external ray for angle ",string(37/63)),
 points(map(realpart,zz4),map(imagpart,zz4)),
 key = concat("external ray for angle ",string(44/63)),
 points(map(realpart,zz5),map(imagpart,zz5)),
 key = concat("external ray for angle ",string(50/63)),
 points(map(realpart,zz6),map(imagpart,zz6)),
 points_joined =false,
 color         = blue,
 point_size    =  0.9,
 key = "repelling fixed point z= beta",
 points(realpart(beta),imagpart(beta)),
 color         = yellow,
 key = "repelling fixed point z= alfa",
 points(realpart(alfa),imagpart(alfa)),
 color         = green,
 key = "period 2 parabolic cycle",
 points(map(realpart,zp),map(imagpart,zp))
);

File history

Click on a date/time to view the file as it appeared at that time.

(Latest | Earliest) View (newer 50) (older 50) (20 | 50 | 100 | 250 | 500)

	Date/Time	Thumbnail	Dimensions	User	Comment
current	18:32, 3 June 2009		1,000×1,000 (22 KB)	Adam majewski (Talk | contribs)	(Julia set and external rays landing on parabolic period 2 cycle)

(Latest | Earliest) View (newer 50) (older 50) (20 | 50 | 100 | 250 | 500)

• Edit this file using an external application (See the setup instructions for more information)

File links