Fractals/Iterations in the complex plane/Julia set
From Wikibooks, the open-content textbooks collection
This book shows how to code different algorithms for drawing sets in dynamical plane : Julia, Filled-in Julia or Fatou sets for complex quadratic polynomial. It is divided in 2 parts :
- description of various algorithms
- descriptions of technics for visualisation of various sets in dynamic plane
- Julia set
- Fatou set
- basin of attraction of infinity ( open set)
- basin of attraction of finite attractor
[edit] Algorithms
[edit] Koenigs coordinate
Koenigs[1] coordinate are used in the basin of attraction of finite attracting (not superattracting) point (cycle),[2][3]
[edit] Escape time
[edit] Basin of attraction to infinity = exterior of filled-in Julia set
Here one computes forward iterations of a complex point Zn:

Here is function which computes the last iteration, that is the first iteration that leaves a circle around the origin with a given escape radius ER for the iteration of the complex quadratic polynomial above. It is a iteration ( integer) for which (abs(Z)>ER).
C version ( here ER2=ER*ER) using double floating point numbers ( without complex type numbers) :
int GiveLastIteration(double Zx, double Zy, double Cx, double Cy, int IterationMax, int ER2) { double Zx2, Zy2; /* Zx2=Zx*Zx; Zy2=Zy*Zy */ int i=0; Zx2=Zx*Zx; Zy2=Zy*Zy; while (i<IterationMax && (Zx2+Zy2<ER2) ) /* ER2=ER*ER */ { Zy=2*Zx*Zy + Cy; Zx=Zx2-Zy2 +Cx; Zx2=Zx*Zx; Zy2=Zy*Zy; i+=1; } return i; }
C++ versions:
int GiveLastIteration(complex C,complex Z , int imax, int ER) { int i; // iteration number for(i=0;i<=imax-1;i++) // forward iteration { Z=Z*Z+C; // overloading of operators if(abs(Z)>ER) break; } return i; }
#include <complex> // C++ complex library // bailout2 = bailout * bailout // this function is based on function esctime from mndlbrot.cpp // from program mandel ver. 5.3 by Wolf Jung // http://www.mndynamics.com/indexp.html int escape_time(complex<double> Z, complex<double> C , int iter_max, double bailout2) { // z= x+ y*i z0=0 long double x =Z.real(), y =Z.imag(), u , v ; int iter; // iteration for ( iter = 0; iter <= iter_max-1; iter++) { u = x*x; v = y*y; if ( u + v <= bailout2 ) { y = 2 * x * y + C.imag(); x = u - v + C.real(); } // if else break; } // for return iter; } // escape_time
Lisp version
((DEFUN GIVELASTITERATION (Z_0 _C IMAX ESCAPE_RADIUS) (SETQ Z Z_0) (SETQ I 0) (LOOP WHILE (AND (< I IMAX) (< (ABS Z) ESCAPE_RADIUS)) DO (INCF I) (SETQ Z (+ (* Z Z) _C))) I)
Maxima version :
GiveLastIteration(z,c):=
block([i:0],
while abs(z)<ER and i<iMax
do (z:z*z + c,i:i+1),
i)$
Delphi version ( using user defined complex type, cabs and f functions )
function GiveLastIteration(z,c:Complex;ER:real;iMax:integer):integer; var i:integer; begin i:=0; while (cabs(z)<ER) and (i<iMax) do begin z:= f(z,c); inc(i); end; result := i; end;
where :
type complex = record x, y: real; end; function cabs(z:complex):real; begin cabs:=sqrt(z.x*z.x+z.y*z.y) end; function f(z,c:complex):complex; // complex quadratic polynomial var tmp:complex; begin tmp.x := (z.x*z.x) - (z.y*z.y) + c.x; tmp.y := 2*z.x*z.y + c.y ; result := tmp; end;
[edit] Boolean Escape time
Algorithm: for every point z of dynamical plane (z-plane) compute iteration number ( last iteration) for which magnitude of z is greater then escape radius. If last_iteration=max_iteration then point is in filled-in Julia set, else it is in its complement(attractive basin of infinity ). Here one has 2 options, so it is named boolean algorithm.
if (LastIteration==IterationMax) then color=BLACK; /* bounded orbits = Filled-in Julia set */ else color=WHITE; /* unbounded orbits = exterior of Filled-in Julia set */
In theory this method is for drawing Filled-in Julia set and its complement ( exterior), but when c is Misiurewicz point ( Filled-in Julia set has no interior) this method draws nothing. For example for c=i . It means that it is good for drawing interior of Filled-in Julia set .
[edit] Integer escape time = Level Sets of the Basin of Attraction of Infinity = Level Sets Method= LSM/J
Escape time measures time of escaping to infinity ( infinity is superattracting point for polynomials). Time is measured in steps ( iterations = i) needed to escape from circle of given radius ( ER= Escape Radius).
One can see few things:
- this is discontinuous function
- i is iMax for z in Filled-in Julia set
- i=0 for x0>ER
- this is nonlinear function
Level sets here are sets of points with the same escape time. Here is algorithm of choosing color in black & white version.
if (LastIteration==IterationMax)
then color=BLACK; /* bounded orbits = Filled-in Julia set */
else /* unbounded orbits = exterior of Filled-in Julia set */
if ((LastIteration%2)==0) /* odd number */
then color=BLACK;
else color=WHITE;
[edit] Normalized iteration count (real escape time)
Math formula :
Maxima version :
GiveNormalizedIteration(z,c,E_R,i_Max):= /* */ block( [i:0,r], while abs(z)<E_R and i<i_Max do (z:z*z + c,i:i+1), r:i-log(log(cabs(z))/log(2)), return(float(r)) )$
In Maxima log(x) is a natural (base e) logarithm of x. To compute log2 use :
log2(x) := log(x) / log(2);
[edit] Level Curves of escape time Method = eLCM/J
These curves are boundaries of Level Sets of escape time ( eLSM/J ). They can be drawn using 2 methods:
- edge detection of Level sets. Algorithm is based on paper by M. Romera et al.[5]
- drawing lemniscates = curves
, see explanation and source code - method described by Harold V. McIntosh[6]
/* Maxima code : draws lemniscates of Julia set */ c: 1*%i; ER:2; z:x+y*%i; f[n](z) := if n=0 then z else (f[n-1](z)^2 + c); load(implicit_plot); /* package by Andrej Vodopivec */ ip_grid:[100,100]; ip_grid_in:[15,15]; implicit_plot(makelist(abs(ev(f[n](z)))=ER,n,1,4),[x,-2.5,2.5],[y,-2.5,2.5]);
[edit] Basin of attraction of finite attractor = interior of filled-in Julia set
- How to find periodic attractor ?
- How many iterations is needed to reach attractor ?
[edit] Components of Interior of Filled Julia set ( Fatou set)
- use limited color ( palette = list of numbered colors)
- find period of attracting cycle
- find one point of attracting cycle
- compute number of iteration after when point reaches the attractor
- color of component=iteration % period[7]
[edit] Decomposition of level sets
[edit] Binary decomposition
Binary decomposition of basin of attraction of infinity (algorithm in pseudocode) :
if (LastIteration==IterationMax)
then color=BLACK; /* bounded orbits = Filled-in Julia set */
else /* unbounded orbits = exterior of Filled-in Julia set */
if (Zy>0) /* Zy=Re(Z) */
then color=BLACK;
else color=WHITE;
[edit] Checking normality
[edit] Checking equicontinuity by Michael Becker
"Iteration is equicontinuous on the Fatou set and not on the Julia set". (Wolf Jung) [8][9]
Michael Becker compares the distance of two close points under iteration on Riemann sphere.[10]
This method can be used to draw not only Julia sets for polynomials ( where infinity is allways superattracting fixed point) but it can be also applied to other functions ( maps), for which infinity is not an attracting fixed point.[11]
[edit] using Marty's criterion by Wolf Jung
Wolf Jung is using "an alternative method of checking normality, which is based on Marty's criterion: |f'| / (1 + |f|^2) must be bounded for all iterates." It is implemented in mndlbrot::marty function ( see src code of program Mandel ver 5.3 ). It uses one point of dynamic plane.
[edit] Inverse Iteration Method (IIM/J) : Julia set
In escape time one computes forward iteration of point z.
In IIM/J one computes:
- repelling fixed point of complex quadratic polynomial

- preimages of
by inverse iterations

Because square root is multivalued function then each
has two preimages
. Thus inverse iteration creates binary tree.
Here is Maxima code for simple IIM:
finverseplus(z,c):=sqrt(z-c); finverseminus(z,c):=-sqrt(z-c); c:%i; /* define c value */ iMax:5000; /* maximal number of reversed iterations */ fixed:float(rectform(solve([z*z+c=z],[z]))); /* compute fixed points of f(z,c) : z=f(z,c) */ if (abs(2*rhs(fixed[1]))<1) /* Find which is repelling */ then block(beta:rhs(fixed[1]),alfa:rhs(fixed[2])) else block(alfa:rhs(fixed[1]),beta:rhs(fixed[2])); z:beta; /* choose repeller as a starting point */ /* make 2 list of points and copy beta to to lists */ xx:makelist (realpart(beta), i, 1, 1); /* list of re(z) */ yy:makelist (imagpart(beta), i, 1, 1); /* list of im(z) */ for i:1 thru iMax step 1 do /* reversed iteration of beta */ block (if random(1.0)>0.5 then z:finverseplus(z,c) else z:finverseminus(z,c), xx:cons(realpart(z),xx), /* save values to draw it later */ yy:cons(imagpart(z),yy) ); plot2d([discrete,xx,yy],[style,[points,1,0,1]]); /* draw reversed orbit of beta */
Compare it with C code.
[edit] Variants of IIM : MIIM/J
- random choose one of two roots ( up to chosen level )= random walk through the tree.
- both roots with the same probability
- more often one then other root
- draw some rare paths in binary tree = MIIM
- draw all roots ( up to chosen level ). There are some posible way to go through all the tree. It is implemented in Fractint[12].
[edit] Complex potential - Boettcher coordinate
Exterior or complement of filled-in Julia set is :
- a basin of attraction of infinity ( superattracting fixed point)
- one of components of the Fatou set

It can be analysed using
- escape time (simple but gives only radial values = escape time ) LSM/J,
- distance estimation ( more advanced, continuus, but gives only radial values = distance ) DEM/J
- Boettcher coordinate or complex potential ( the best , gives :
- radial values ( real potential ) LogPhi = CPM/J
- angular values ( external angle ) ArgPhi
So both values can be used to color with 2D gradient.
First compute Boettcher coordinate
using Wolf Jung formula [13]

[edit] LogPhi - Douady-Hubbard potential - real potential - radial component of complex potential
Note that potential inside Kc is zero so (in pseudocode ) :
if (LastIteration==IterationMax) then potential=0 /* inside Filled-in Julia set */ else potential= GiveLogPhi(z0,c,ER,nMax); /* outside */
It also removes potential error for log(0).
Math (full) version : [14]

Maxima (full) version :
GiveLogPhi(z0,c,ER,nMax):= block( [z:z0, logphi:log(cabs(z)), fac:1/2, n:0], while n<nMax and abs(z)<ER do (z:z*z+c, logphi:logphi+fac*log(cabs(1+c/(z*z))), n:n+1 ), return(float(logphi)) )$
Math simplified formula :

Maxima function :
GiveSLogPhi(z0,c,e_r,i_max):=
block(
[z:z0,
logphi,
fac:1/2,
i:0
],
while i<i_max and cabs(z)<e_r do
(z:z*z+c,
fac:fac/2,
i:i+1
),
logphi:fac*log(cabs(z)),
return(float(logphi))
)$
If you dont check if orbit is not bounded ( escapes, bailout test) then use this Maxima function :
GiveSLogPhi(z0,c,e_r,i_max):=
block(
[z:z0, logphi, fac:1/2, i:0],
while i<i_max and cabs(z)<e_r do
(z:z*z+c,
fac:fac/2,
i:i+1 ),
if i=i_max
then logphi:0
else logphi:fac*log(cabs(z)),
float(logphi)
)$
C version :
double jlogphi(double zx0, double zy0, double cx, double cy)
/* this function is based on function by W Jung */
{
int j;
double
zx=zx0,
zy=zy0,
s = 0.5,
zx2=zx*zx,
zy2=zy*zy,
t;
for (j = 1; j < 400; j++)
{ s *= 0.5;
zy = 2 * zx * zy + cy;
zx = zx2 - zy2 + cx;
zx2 = zx*zx;
zy2 = zy*zy;
t = fabs(zx2 + zy2); // abs(z)
if ( t > 1e24) break;
}
return s*log2(t); // log(zn)* 2^(-n)
}//jlogphi
[edit] CPM/J
[edit] Level Sets of potential = pLSM/J
Here is Delphi function which gives level of potential :
Function GiveLevelOfPotential(potential:extended):integer;
var r:extended;
begin
r:= log2(abs(potential));
result:=ceil(r);
end;
[edit] Level Curves of potential = equipotential lines = pLCM/J
[edit] ArgPhi - External angle - angular component of complex potential
One can start with binary decomposition of basin of attraction of infinity.
The second step can be using 
Computing external angle
External angle (argument) is argument of Boettcher coordinate 

Because Boettcher coordinate is a product of complex numbers so argument of product is :

[edit] Constructing the spine of filled Julia set
Algorithm for constructiong the spine is described by A. Douady[15]
- join
and
, - (to do )
[edit] Drawing dynamic external ray
[edit] backwards iteration
This method have been used by several people and proved by Thierry Bousch. [16]
Code in c++ by Wolf Jung can be found in procedure QmnPlane::backray() in file qmnplane.cpp ( see source code of program mandel version 5.3 ). [17]
- Periodic ray ( simplest case )
It will be explained by an example :
First choose external angle
(in turns). External angle for periodic ray is a rational number.
Compute period of external angle under doubling map.
Because "1/3 doubled gives 2/3 and 2/3 doubled gives 4/3, which is congruent to 1/3" [18]
or
so external angle
has period 2 under doubling map.
Start with 2 points near infinity (in conjugate plane): on ray 1/3
and on ray 2/3
, here
so one can swith to dynamical plane ( Boettcher conjugation )
Backward iteration (with proper chose from two possibilities)[20] of point on ray 1/3 goes to ray 2/3, back to 1/3 and so on.
Backward iteration (with proper chose from two possibilities) of point on ray 2/3 goes to ray 1/3, back to 2/3 and so on.
One has to divide points of backward orbit between 2 rays and can draw this 2 rays.
- Preperiodic ray
[edit] Drawing dynamic external ray using inverse Boettcher map by Curtis McMullen
This method is based on C program by Curtis McMullen[21] and its Pascal version by Matjaz Erat[22]
It consist of 3 big steps :
- compute some w-points of external ray of circle for angle
and various radii (rasterisation)
where 
- map w-points to z-point using inverse Boettcher map

- draw z-points ( and conect them using lines )
First and last steps are easy, but second is not so needs more explanation.
[edit] Rasterisation
For given external ray in
plane each point of ray has :
- constant value
( external angle in turns ) - variable radius

so
points of ray are parametrised by radius
and can be computed using exponential form of complex numbers :
One can go along ray using linear scale :
t:1/3; /* example value */ R_Max:4; R_Min:1.1; for R:R_Max step -0.5 thru R_Min do w:R*exp(2*%pi*%i*t); /* Maxima allows non-integer values in for statement */
It gives some w points with equal distance between them.
Another method is to use nonlinera scale.
To do it we introduce floating point exponent
such that :
and
To compute some w points of external ray in
plane for angle
use such Maxima code :
t:1/3; /* external angle in turns */ /* range for computing R ; as r tends to 0 R tends to 1 */ rMax:2; /* so Rmax=2^2=4 / rMin:0.1; /* rMin > 0 */ caution:0.93; /* positive number < 1 ; r:r*caution gives smaller r */ r:rMax; unless r<rMin do ( r:r*caution, /* new smaller r */ R:2^r, /* new smaller R */ w:R*exp(2*%pi*%i*t) /* new point w in f0 plane */ );
In this method distance between points is not equal but inversely proportional to distance to boundary of filled Julia set.
It is good because here ray has greater curvature so curve will be more smooth.
[edit] Mapping
Mapping points from
-plane to
-plane consist of 4 minor steps :
- forward iteration in
plane
until
is near infinity
- switching plane ( from
to
)
( because here, near infinity :
)
- backward iteration in
plane the same (
) number of iterations - last point
is on our external ray
1,2 and 4 minor steps are easy. Third is not.
Backward iteration uses square root of complex number. It is 2-valued functions so backward iteration gives binary tree.
One can't choose good path in such tree without extre informations. To solve it we will use 2 things :
- equicontinuity of basin of attraction of infinity
- conjugacy between
and
planes
[edit] =Equicontinuity of basin of attraction of infinity=
Basin of attraction of infinity ( complement of filled-in Julia set) contains all points which tends to infinity under forward iteration.

Infinity is superattracting fixed point and orbits of all points have similar behaviour. In other words orbits of 2 points are assumed to stay close if they are close at the beginning.
It is equicontinuity ( compare with normality).
In
plane one can use forward orbit of previous point of ray for computing backward orbit of next point.
[edit] Detailed version of algorithm
- compute first point of ray (start near infinity ang go toward Julia set )
where 
here one can easily switch planes :
It is our first z-point of ray.
- compute next z-point of ray
- compute next w-point of ray for

- compute forward iteration of 2 points : previous z-point and actual w-point. Save z-orbit and last w-point
- switch planes and use last w-point as a starting point :

- backward iteration of new
toward new
using forward orbit of previous z point
is our next z point of our ray
- compute next w-point of ray for
- and so on ( next points ) until

[edit] Maxima CAS src code
/* gives a list of z-points of external ray for angle t in turns and coefficient c */ 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) )$
[edit] BSM/J
[edit] DEM/J
[edit] Internal distance estimation
[edit] External distance estimation
" For distance estimate it has been proved that the computed value differs from the true distance at most by a factor of 4. " (Wolf Jung)
Math formula :

where :
is first derivative with respect to c.
This derivative can be found by iteration starting with
and then :

Algorithm for computing distance from point z ( in exterior ) to nearest point on the boundary of Julia set. Here is pseudocode:
if (LastIteration==IterationMax)
then { /* interior of Julia set, distance = 0 , color black */ }
else /* exterior of Filled-in Julia set */
{ double distance=give_distance(Z0,C,IterationMax);
if (distance<distanceMax)
then { /* Julia set : color = white */ }
else { /* exterior of Julia set : color = black */}
}
DistanceMax is smaller then pixel size. One can start with DistanceMax=pixel_size and if the result is not good change DistanceMax=pixel_size/n where n is number greater then 1.
Examples of code :
- in cpp see mndlbrot::dist from mndlbrot.cpp in src code of program mandel ver 5.3 by Wolf Jung.
- Maxima function
GiveExtDistance(z0,c,e_r,i_max):= /* needs z in exterior of Kc */ block( [z:z0, dz:1, cabsz:cabs(z), cabsdz:1, /* overflow limit */ i:0], while cabsdz < 10000000 and i<i_max do ( dz:2*z*dz, z:z*z + c, cabsdz:cabs(dz), i:i+1 ), cabsz:cabs(z), return(2*cabsz*log(cabsz)/cabsdz) )$
[edit] Color
- Visualising Julia sets by Georg-Johann
- Combined Methods of Depicting Julia Sets and Parameter Planes by Chris King
- On Fractal Coloring Techniques by Jussi Harkonen Master's Thesis, Department of Mathematics, Åbo Akademi University, Turku, 2007, 61 pages. The thesis was carried out under the supervision of Professor Goran Hognas
- Technical Info - Colorizing by Michael Condron
- Colors/Color_gradient
[edit] Sets
[edit] Julia sets
"Most programs for computing Julia sets work well when the underlying dynamics is hyperbolic but experience an exponential slowdown in the parabolic case." ( Mark Braverman )[23]
- IIM/J
- DEM/J
[edit] Fatou set
[edit] Periodic points
[edit] Finding equations for periodic points
[edit] Standard equations
Equation for periodic point for period
is :

so in Maxima :
for i:1 thru 4 step 1 do disp(i,expand(Fn(i,z,c)=0)); 1 z^2-z+c=0 2 z^4+2*c*z^2-z+c^2+c=0 3 z^8+4*c*z^6+6*c^2*z^4+2*c*z^4+4*c^3*z^2+4*c^2*z^2-z+c^4+2*c^3+c^2+c=0 4 z^16+8*c*z^14+28*c^2*z^12+4*c*z^12+56*c^3*z^10+24*c^2*z^10+70*c^4*z^8+60*c^3*z^8+6*c^2*z^8+2*c*z^8+56*c^5*z^6+80*c^4*z^6+ 24*c^3*z^6+8*c^2*z^6+28*c^6*z^4+60*c^5*z^4+36*c^4*z^4+16*c^3*z^4+4*c^2*z^4+8*c^7*z^2+24*c^6*z^2+24*c^5*z^2+16*c^4*z^2+ 8*c^3*z^2- z+c^8+4*c^7+6*c^6+6*c^5+5*c^4+2*c^3+c^2+c=0
Degree of standard equations is 
These equations give periodic points for period p and its divisors. It is because these polynomials are product of lower period polynomials
:
for i:1 thru 4 step 1 do disp(i,factor(Fz(i,z,c))); 1 z^2-z+c 2 (z^2-z+c)*(z^2+z+c+1) 3 (z^2-z+c)*(z^6+z^5+3*c*z^4+z^4+2*c*z^3+z^3+3*c^2*z^2+3*c*z^2+z^2+c^2*z+2*c*z+z+c^3+2*c^2+c+1) 4 (z^2-z+c)*(z^2+z+c+1)(z^12+6*c*z^10+z^9+15*c^2*z^8+3*c*z^8+4*c*z^7+20*c^3*z^6+12*c^2*z^6+z^6+6*c^2*z^5+2*c*z^5+15*c^4*z^4+18*c^3*z^4+3* c^2*z^4+4*c*z^4+4*c^3*z^3+4*c^2*z^3+z^3+6*c^5*z^2+12*c^4*z^2+6*c^3*z^2+5*c^2*z^2+c*z^2+c^4*z+2*c^3*z+c^2*z+2*c*z+c^6+ 3*c^5+3*c^4+3*c^3+2*c^2+1)
so





Standard equations can be reduced to equations
giving periodic points for period p without it's divisors.
[edit] Reduced equations

so in Maxima gives :
for i:1 thru 5 step 1 do disp(i,expand(G[i,z,c]=0)); 1 z^2-z+c=0 2 z^2+z+c+1=0 3 z^6+z^5+3*c*z^4+z^4+2*c*z^3+z^3+3*c^2*z^2+3*c*z^2+z^2+c^2*z+2*c*z+z+c^3+2*c^2+c+1=0 4 z^12+6*c*z^10+z^9+15*c^2*z^8+3*c*z^8+4*c*z^7+20*c^3*z^6+12*c^2*z^6+z^6+6*c^2*z^5+2*c*z^5+15*c^4*z^4+18*c^3*z^4+ 3*c^2*z^4+4*c*z^4+4*c^3*z^3+4*c^2*z^3+z^3+6*c^5*z^2+12*c^4*z^2+ 6*c^3*z^2+5*c^2*z^2+c*z^2+c^4*z+2*c^3*z+c^2*z+2*c*z+c^6+3*c^5+3*c^4+3*c^3+2*c^2+1=0
[edit] Computing periodic points in Maxima
Define c value
(%i1) c:%i; (%o1) %i
[edit] Fixed points ( period = 1 )
Compute fixed points
(%i2) p:float(rectform(solve([z*z+c=z],[z]))); (%o2) [z=0.62481053384383*%i-0.30024259022012,z=1.30024259022012-0.62481053384383*%i]
Find which is repelling :
abs(2*rhs(p[1]));
Show that sum of fixed points is 1
(%i10) p:solve([z*z+c=z], [z]); (%o10) [z=-(sqrt(1-4*%i)-1)/2,z=(sqrt(1-4*%i)+1)/2] (%i14) s:radcan(rhs(p[1]+p[2])); (%o14) 1
Draw points :
(%i15) xx:makelist (realpart(rhs(p[i])), i, 1, length(p)); (%o15) [-0.30024259022012,1.30024259022012] (%i16) yy:makelist (imagpart(rhs(p[i])), i, 1, length(p)); (%o16) [0.62481053384383,-0.62481053384383] (%i18) plot2d([discrete,xx,yy], [style, points]);
One can add point z=1/2 to the lists:
(%i31) xx:cons(1/2,xx); (%o31) [1/2,-0.30024259022012,1.30024259022012] (%i34) yy:cons(0,yy); (%o34) [0,0.62481053384383,-0.62481053384383] (%i35) plot2d([discrete,xx,yy], [style, points]);
[edit] period 2 points
(%i2) solve([(z*z+c)^2+c=z], [z]); (%o2) [z=-(sqrt(-4*c-3)+1)/2,z=(sqrt(-4*c-3)-1)/2,z=-(sqrt(1-4*c)-1)/2,z=(sqrt(1-4*c)+1)/2]
Show that z1+z2 = -1
(%i4) s:radcan(rhs(p[1]+p[2])); (%o4) -1
Show that z2+z3=1
(%i6) s:radcan(rhs(p[3]+p[4])); (%o6) 1
[edit] Finding period of orbit
[edit] More tutorials and code
- in Java see Evgeny Demidov
- in C see :
- in C++ see Wolf Jung page,
- in Gnuplot see Tutorial by T.Kawano
- in Lisp for Maxima see Dynamics by Jaime E. Villate
- in Mathemathica see :
[edit] References
- ↑ Gabriel Koenigs biographie at The MacTutor History of Mathematics archive
- ↑ G. Koenigs, Recherches sur les intégrales de certaines équations fonctionnelles, Annales École Normale Supérieure, 1(3) (1884), Supplément, 3-41.
- ↑ Inigo Quilez images and tutuorial
- ↑ Program Mandel by Wolf Jung on GNU General Public License
- ↑ Drawing the Mandelbrot set by the method of escape lines. M. Romera et al.
- ↑ Julia Curves, Mandelbrot Set, Harold V. McIntosh.
- ↑ The fixed points and periodic orbits by E Demidov
- ↑ Alan F. Beardon, S. Axler, F.W. Gehring, K.A. Ribet : Iteration of Rational Functions: Complex Analytic Dynamical Systems. Springer, 2000; ISBN 0387951512, 9780387951515; page 49
- ↑ Joseph H. Silverman : The arithmetic of dynamical systems. Springer, 2007. ISBN 0387699031, 9780387699035; page 22
- ↑ Visualising Julia sets by Georg-Johann
- ↑ Julia sets by Michael Becker. See the metric d(z,w)
- ↑ Fractint documentation - Inverse Julias
- ↑ How to draw external rays by Wolf Jung
- ↑ The Beauty of Fractals, page 65
- ↑ A. Douady, “Algorithms for computing angles in the Mandelbrot set,” in Chaotic Dynamics and Fractals, M. Barnsley and S. G. Demko, Eds., vol. 2 of Notes and Reports in Mathematics in Science and Engineering, pp. 155–168, Academic Press, Atlanta, Ga, USA, 1986.
- ↑ Thierry Bousch : De combien tournent les rayons externes? Manuscrit non publié, 1995
- ↑ Program Mandel by Wolf Jung
- ↑ Explanation by Wolf Jung
- ↑ Modular arithmetic in wikipedia
- ↑ Square root of complex number gives 2 values so one has to choose only one. For details see Wolf Jung page
- ↑ c program by Curtis McMullen (quad.c in Julia.tar.gz)
- ↑ Quadratische Polynome by Matjaz Erat
- ↑ Mark Braverman : On efficient computation of parabolic Julia sets










