# File:JuliaRay 2 5.png

Julia set and external rays landing on parabolic period 2 cycle

# Algorithm

Program draws (to png file) :

• dynamic rays $\mathcal{R}_t = {z : arg_e(z)=t } \,$ using $z = \Psi_n(w) = f_c^{-n}(w^{2^n})\,$
• Julia set ( backward orbit of repeller = IIM/J )

## How to compute $c\,$ coefficient ?

We want point c for which period 2 cycle z0,z1 is parabolic ( indifferent ). It is when absolute value of multiplier at periodic point $z_0\,$ is equal to one :

First compute point on unit circle with internal angle $\frac{1}{3}\,$ ( any other rational angle will also be good) :

$w = e^{\frac{2\pi i}{3} }\,$

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 $\gamma_2\,$:

$\gamma_2(w) = \frac{w}{4} - 1\,$

so

$c = \frac{e^{\frac{2\pi i}{3} }}{4} - 1 = -1.125 +0.21650635094611*i\,$

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 */
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)

The following page links to this file: