Fractals/Iterations in the complex plane/Mandelbrot set
From Wikibooks, the open-content textbooks collection
This book shows how to code different algorithms for drawing parameter plane (Mandelbrot set ) for complex quadratic polynomial.
[edit] Exterior of Mandelbrot set
[edit] Target set
Target set
is an arbitrary set on dynamical plane containing infinity and not containing points of Filled-in Fatou sets.
For escape time algorithms target set determines the shape of level sets and curves. It does not do it for other methods.
- Exterior of circle
This is typical target set. It is exterior of circle with center at origin
and radius =ER :

Radius is named escape radius ( ER ) or bailout value. Radius should be greater than 2.
- Exterior of square
Here target set is exterior of square of side length
centered at origin

[edit] Bailout test
[edit] Escape time
[edit] Boolean escape time
Here complex plane consists of 2 sets : Mandelbrot set
and its complement
:

This is complete code of C program.
- It makes a ppm file which consists an image. To see the file (image) use external application ( graphic viewer).
- Program consists of 3 loops:
- iY and iX, which are used to scan rectangle area of parameter plane
- iterations.
For each point of screen (iX,iY) it's complex value is computed c=cx+cy*i.
For each point c is computed iterations of critical point 
...
/* c program: -------------------------------- 1. draws Mandelbrot set for Fc(z)=z*z +c using Mandelbrot algorithm ( boolean escape time ) ------------------------------- 2. technique of creating ppm file is based on the code of Claudio Rocchini http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg create 24 bit color graphic file , portable pixmap file = PPM see http://en.wikipedia.org/wiki/Portable_pixmap to see the file use external application ( graphic viewer) */ #include <stdio.h> int main() { /* screen ( integer) coordinate */ int iX,iY; const int iXmax = 800; const int iYmax = 800; /* world ( double) coordinate = parameter plane*/ double Cx,Cy; const double CxMin=-2.5; const double CxMax=1.5; const double CyMin=-2.0; const double CyMax=2.0; /* */ double PixelWidth=(CxMax-CxMin)/iXmax; double PixelHeight=(CyMax-CyMin)/iYmax; /* color component ( R or G or B) is coded from 0 to 255 */ /* it is 24 bit color RGB file */ const int MaxColorComponentValue=255; FILE * fp; char *filename="new1.ppm"; char *comment="# ";/* comment should start with # */ static unsigned char color[3]; /* Z=Zx+Zy*i ; Z0 = 0 */ double Zx, Zy; double Zx2, Zy2; /* Zx2=Zx*Zx; Zy2=Zy*Zy */ /* */ int Iteration; const int IterationMax=200; /* bail-out value , radius of circle ; */ const double EscapeRadius=2; double ER2=EscapeRadius*EscapeRadius; /*create new file,give it a name and open it in binary mode */ fp= fopen(filename,"wb"); /* b - binary mode */ /*write ASCII header to the file*/ fprintf(fp,"P6\n %s\n %d\n %d\n %d\n",comment,iXmax,iYmax,MaxColorComponentValue); /* compute and write image data bytes to the file*/ for(iY=0;iY<iYmax;iY++) { Cy=CyMin + iY*PixelHeight; if (fabs(Cy)< PixelHeight/2) Cy=0.0; /* Main antenna */ for(iX=0;iX<iXmax;iX++) { Cx=CxMin + iX*PixelWidth; /* initial value of orbit = critical point Z= 0 */ Zx=0.0; Zy=0.0; Zx2=Zx*Zx; Zy2=Zy*Zy; /* */ for (Iteration=0;Iteration<IterationMax && ((Zx2+Zy2)<ER2);Iteration++) { Zy=2*Zx*Zy + Cy; Zx=Zx2-Zy2 +Cx; Zx2=Zx*Zx; Zy2=Zy*Zy; }; /* compute pixel color (24 bit = 3 bajts) */ if (Iteration==IterationMax) { /* interior of Mandelbrot set = black */ color[0]=0; color[1]=0; color[2]=0; } else { /* exterior of Mandelbrot set = white */ color[0]=255; /* Red*/ color[1]=255; /* Green */ color[2]=255;/* Blue */ }; /*write color to the file*/ fwrite(color,1,3,fp); } } fclose(fp); return 0; }
[edit] Integer escape time = LSM/M
This is also called Level Set Method ( LSM )

Diffrence between Mandelbrot algorithm and LSM/M is in only in part instruction , which computes pixel color of exterior of Mandelbrot set. In LSM/M is :
if (Iteration==IterationMax) { /* interior of Mandelbrot set = black */ color[0]=0; color[1]=0; color[2]=0; } /* exterior of Mandelbrot set = LSM */ else if ((Iteration%2)==0) { /* even number = black */ color[0]=0; /* Red*/ color[1]=0; /* Green */ color[2]=0;/* Blue */ } else {/* odd number = white */ color[0]=255; /* Red*/ color[1]=255; /* Green */ color[2]=255;/* Blue */ };
[edit] Optimisation
Here is C++ function which can be used to draw LSM/M :
int iterate_mandel(complex C , int imax, int bailout) { int i; std::complex Z(0,0); // initial value for iteration Z0 for(i=0;i<=imax-1;i++) { Z=Z*Z+C; // overloading of operators if(abs(Z)>bailout)break; } return i; }
I think that it can't be coded simpler (it looks better then pseudocode), but it can be coded in other thay which can be executed faster .
Here is faster C function which can be used instead of above C++ function:
int GiveEscapeTime(double C_x, double C_y, int iMax, double _ER2) { int i; double Zx, Zy; double Zx2, Zy2; /* Zx2=Zx*Zx; Zy2=Zy*Zy */ Zx=0.0; /* initial value of orbit = critical point Z= 0 */ Zy=0.0; Zx2=Zx*Zx; Zy2=Zy*Zy; for (i=0;i<iMax && ((Zx2+Zy2)<_ER2);i++) { Zy=2*Zx*Zy + C_y; Zx=Zx2-Zy2 +C_x; Zx2=Zx*Zx; Zy2=Zy*Zy; }; return i; }
Above function doesn't use explicit definition of complex number.
[edit] Level Curves of escape time Method = LCM/M
Lemniscates are boundaries of Level Sets of escape time ( LSM/M ). They can be drawn using 2 methods:
- edge detection of Level sets. Algorithm is based on paper by M. Romera et al.[1] This method is fast and allows looking for high iterations.
- boundary trace[2]
- drawing curves
, see explanation and source code. This method is very complicated for iterations>5.
[edit] Decomposition of exterior of Mandelbrot set
Decomposition is modification of escape time algorithm.
The target set is divided into parts (2 or more). Very large escape radius is used, for example ER = 12.
[edit] Binary decomposition of LSM/M
Here target set
on dynamic plane is divided into 2 parts (binary decomposition = 2-decomposition ):
- upper half ( white)

- lower half (black)

Division of target set induces decomposition of level sets
into
parts:
which is colored white,
which is colored black.
External rays of angles (measured in turns):

can be seen.
Difference between binary decomposition algorithm and Mandel or LSM/M is in only in part of instruction , which computes pixel color of exterior of Mandelbrot set. In binary decomposition is :
if (Iteration==IterationMax) { /* interior of Mandelbrot set = black */ color[0]=0; color[1]=0; color[2]=0; } /* exterior of Mandelbrot set = LSM */ else if (Zy>0) { color[0]=0; /* Red*/ color[1]=0; /* Green */ color[2]=0;/* Blue */ } else { color[0]=255; /* Red*/ color[1]=255; /* Green */ color[2]=255;/* Blue */ };
Point c is plotting white or black if imaginary value of last iteration ( Zy) is positive or negative.[3]
[edit] Distance estimation DEM/M
Estimates distance from point
( in the exterior of M ) to nearest point in M.
Math formula : [4]
where :
[edit] Complex potential
Complex potential is a complex number, so it has a real (real potential) and an imaginary part(external angle). One can take its curl, divergence or it's absolute value. So on one image one can use more then one variable to color image.[5]
[edit] Real potential = CPM/M
One can use potential to:
- smooth ( continuus) coloring
- discrete coloring ( level sets of potential)
- 3D view
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] External angle and external ( parameter) ray
[edit] Boundary of Mandelbrot set
Methods used to draw boundary of Mandelbrot set :
- DEM/M
- Boundary scanning : detecting the edges after Boolean escape time
- contour integration by Robert Davies[8]
- Jungreis function

[edit] Computing Misiurewicz points of complex quadratic mapping in Maxima
Misiurewicz points are special boundary points.
Define polynomial.
P(n):=if n=0 then 0 else P(n-1)^2+c;
Define a function whose roots are Misiurewicz points, and find them.
M(preperiod,period):=allroots(%i*P(preperiod+period)-%i*P(preperiod));
Examples of use :
(%i6) M(2,1); (%o6) [c=-2.0,c=0.0] (%i7) M(2,2); (%o7) [c=-1.0*%i,c=%i,c=-2.0,c=-1.0,c=0.0]
[edit] Interior of Mandelbrot set - hyperbolic components
[edit] Interior distance estimation
Description of method [9]
[edit] absolute value of the orbit
[edit] bof60
This is the method described in the book "The Beauty of Fractals" on page 60. It is used only for interior points of the Mandelbrot set.
Color of point is proportional to smallest distance of its orbit from origin[10]. Level sets of distance are sets of points with the same distance[11]
[edit] bof61
This is the method described in the book "The Beauty of Fractals" on page 61.
Index (c) is the iteration when point of the orbit was closest to the origin. Since there may be more than one, index(c) is the least such. This algorithms shows borders of domains with the same index(c)[12] [13].
[edit] Period of hyperbolic components
[edit] Multiplier map
[edit] Internal angle
[edit] Internal rays
[edit] Boundaries of hyperbolic components
Methods of drawing boundaries:
- solving boundary equations ( Brown,Stephenson,Jung ; works only up to period 7-8 )
- parametrisation of boundary with Newton method near centers of components [14]. This methods needs centers, so it has some limitations,
- Boundary scanning : detecting the edges after detecting period by checking every pixel. This method is slow but allows zooming. Draws "all" components
- Boundary tracing for given c value. Draws single component.
- Fake Mandelbrot set by Anne M. Burns : draws main cardioid and all its descendants. Do not draw mini Mandelbrot sets. [15]
[edit] solving boundary equations
[edit] System of 2 equations defining boundaries of period
hyperbolic components
- first defines periodic point,
- second defines indifferent orbit.

Because stability index
is equal to radius of point of unit circle
:

so one can change second equation to form [16] :

It gives system of equations :

Above system of 2 equations has 3 variables :
(
is constant and multiplier
is a function of
). One have to remove 1 variable to be able to solve it.
Boundaries are closed curves : cardioids or circles. One can parametrize points of boundaries with angle
( here measured in turns from 0 to 1 ).
After evaluation of
one can put it into above system, and get a system of 2 equations with 2 variables
.
Now it can be solved
For periods:
- 1-3 it can be solved with symbolical methods and give implicit ( boundary equation)
and explicit function (inverse multiplier map) :
- 1-2 it is easy to solve [17]
- 3 it can be solve using "elementary algebra" ( Stephenson )
- >3 it can't be reduced to explicitly function but
- can be reduced to implicit solution ( boundary equation)
and solved numerically - can be solved numerically using Newton method for solving system of nonlinear equations
- can be reduced to implicit solution ( boundary equation)
[edit] Solving system of equation for period 1
Here is Maxima code :
(%i4) p:1; (%o4) 1 (%i5) e1:F(p,z,c)=z; (%o5) z^2+c=z (%i6) e2:m(p)=w; (%o6) 2*z=w (%i8) s:eliminate ([e1,e2], [z]); (%o8) [w^2-2*w+4*c] (%i12) s:solve([s[1]], [c]); (%o12) [c=-(w^2-2*w)/4] (%i13) define (m1(w),rhs(s[1])); (%o13) m1(w):=-(w^2-2*w)/4
Second equation contains only one variable, one can eliminate this variable. Because boundary equation is simple so it is easy to get explicit solution
m1(w):=-(w^2-2*w)/4
[edit] Solving system of equation for period 2
Here is Maxima code using to_poly_solve package by Barton Willis:
(%i4) p:2; (%o4) 2 (%i5) e1:F(p,z,c)=z; (%o5) (z^2+c)^2+c=z (%i6) e2:m(p)=w; (%o6) 4*z*(z^2+c)=w (%i7) e1:F(p,z,c)=z; (%o7) (z^2+c)^2+c=z (%i10) load(topoly_solver); to_poly_solve([e1, e2], [z, c]); (%o10) C:/PROGRA~1/MAXIMA~1.1/share/maxima/5.16.1/share/contrib/topoly_solver.mac (%o11) [[z=sqrt(w)/2,c=-(w-2*sqrt(w))/4],[z=-sqrt(w)/2,c=-(w+2*sqrt(w))/4],[z=(sqrt(1-w)-1)/2,c=(w-4)/4],[z=-(sqrt(1-w)+1)/2,c=(w-4)/4]] (%i12) s:to_poly_solve([e1, e2], [z, c]); (%o12) [[z=sqrt(w)/2,c=-(w-2*sqrt(w))/4],[z=-sqrt(w)/2,c=-(w+2*sqrt(w))/4],[z=(sqrt(1-w)-1)/2,c=(w-4)/4],[z=-(sqrt(1-w)+1)/2,c=(w-4)/4]] (%i14) rhs(s[4][2]); (%o14) (w-4)/4 (%i16) define (m2 (w),rhs(s[4][2])); (%o16) m2(w):=(w-4)/4
explicit solution :
m2(w):=(w-4)/4
[edit] Solving system of equation for period 3
For period 3 ( and higher) previous method give no results (Maxima code) :
(%i14) p:3; e1:z=F(p,z,c); e2:m(p)=w; load(topoly_solver); to_poly_solve([e1, e2], [z, c]); (%o14) 3 (%o15) z=((z^2+c)^2+c)^2+c (%o16) 8*z*(z^2+c)*((z^2+c)^2+c)=w (%i17) (%o17) C:/PROGRA~1/MAXIMA~1.1/share/maxima/5.16.1/share/contrib/topoly_solver.mac `algsys' cannot solve - system too complicated. #0: to_poly_solve(e=[z = ((z^2+c)^2+c)^2+c,8*z*(z^2+c)*((z^2+c)^2+c) = w],vars=[z,c]) -- an error. To debug this try debugmode(true);
I use code by Robert P. Munafo[18] which is based on paper of Wolf Jung.
[edit] Boundary equation
The result of solving above system with respect to
is boundary equation,
where
is boundary polynomial.
It defines exact coordinates of hyperbolic componenets for given period
.
It is implicit equation.
[edit] Solving boundary equations
Solving boundary euations for various angles gives list of boundary points.
[edit] Centers of components
Methods of finding centers :
- "all" centers
- scanning parameter plane for ponts c for which orbit of critical point vanishes [19]
- all centers for given period
- one center
for given period
near given point
- some steps of Newton method for solving

- some steps of Newton method for solving
[edit] System of equations defining centers
Center of hyperbolic component is a point
of parameter plain, for which periodic z-cycle is superattracting. It gives a system of 2 equations defining centers of period
hyperbolic components :
- first defines periodic point
, - second makes
point superattracting.

see definitions for details.
[edit] Solving system of equations
Second equation can be solved when critical point
is in this cycle :

To solve system put
into first equation.
[edit] Equation defining centers
One gets equation :

Centers of components are roots of above equation.
[edit] Reduced equation defining centers
Set of roots of above equation contains centers for period p and its divisors. It is because :

where :
is irreducible divisors of
[22][1]- capital Pi notation of iterated multiplication is used
means here : for all divisors of
( from 1 to p ). See table of divisors.
For example :




So one can find irreducible polynomials using :

Here is Maxima function for computing
:
GiveG[p]:= block( [f:divisors(p), t:1], /* t is temporary variable = product of Gn for (divisors of p) other then p */ f:delete(p,f), /* delete p from list of divisors */ if p=1 then return(Fn(p,0,c)), for i in f do t:t*GiveG[i], g: Fn(p,0,c)/t, return(ratsimp(g)) )$
Here is a table with degree of
and
for periods 1-10 and precision needed to compute roots of these functions.
period ![]() |
![]() |
![]() |
![]() |
|---|---|---|---|
| 1 | 1 | 1 | 16 |
| 2 | 2 | 1 | 16 |
| 3 | 4 | 3 | 16 |
| 4 | 8 | 6 | 16 |
| 5 | 16 | 15 | 16 |
| 6 | 32 | 27 | 16 |
| 7 | 64 | 63 | 32 |
| 8 | 128 | 120 | 64 |
| 9 | 256 | 252 | 128 |
| 10 | 512 | 495 | 300 |
| 11 | 1024 | 1023 | 600 |
Here is a table of greatest coefficients.
period ![]() |
![]() |
![]() |
|---|---|---|
| 1 | 0 | 1 |
| 2 | 0 | 1 |
| 3 | 1 | 2 |
| 4 | 1 | 3 |
| 5 | 3 | 116 |
| 6 | 4 | 5892 |
| 7 | 11 | 17999433372 |
| 8 | 21 | 106250977507865123996 |
| 9 | 44 | 16793767022807396063419059642469257036652437 |
| 10 | 86 | 86283684091087378792197424215901018542592662739248420412325877158692469321558575676264 |
| 11 | 179 | 307954157519416310480198744646044941023074672212201592336666825190665221680585174032224052483643672228833833882969097257874618885560816937172481027939807172551469349507844122611544 |
Precision can be estimated as bigger then size of binary form of greatest coefficient
:
period ![]() |
![]() |
![]() |
![]() |
|---|---|---|---|
| 1 | 1 | 1 | 16 |
| 2 | 1 | 1 | 16 |
| 3 | 3 | 2 | 16 |
| 4 | 6 | 2 | 16 |
| 5 | 15 | 7 | 16 |
| 6 | 27 | 13 | 16 |
| 7 | 63 | 34 | 32 |
| 8 | 120 | 67 | 64 |
| 9 | 252 | 144 | 128 |
| 10 | 495 | 286 | 300 |
| 11 | 1023 | 596 | 600 |
Here is Maxima code for 
period_Max:11; /* ----------------- definitions -------------------------------*/ /* basic funtion = monic and centered complex quadratic polynomial http://en.wikipedia.org/wiki/Complex_quadratic_polynomial */ f(z,c):=z*z+c $ /* iterated function */ fn(n, z, c) := if n=1 then f(z,c) else f(fn(n-1, z, c),c) $ /* roots of Fn are periodic point of fn function */ Fn(n,z,c):=fn(n, z, c)-z $ /* gives irreducible divisors of polynomial Fn[p,z=0,c] */ GiveG[p]:= block( [f:divisors(p),t:1], g, f:delete(p,f), if p=1 then return(Fn(p,0,c)), for i in f do t:t*GiveG[i], g: Fn(p,0,c)/t, return(ratsimp(g)) )$ /* degree of function */ GiveDegree(_function,_var):=hipow(expand(_function),_var); log10(x) := log(x)/log(10); /* ------------------------------*/ file_output_append:true; /* to append new string to file */ grind:true; for period:1 thru period_Max step 1 do ( g[period]:GiveG[period], /* function g */ d[period]:GiveDegree(g[period],c), /* degree of function g */ cf[period]:makelist(coeff(g[period],c,d[period]-i),i,0,d[period]), /* list of coefficients */ cf_max[period]:apply(max, cf[period]), /* max coefficient */ disp(cf_max[period]," ",ceiling(log10(cf_max[period]))), s:concat("period:",string(period)," cf_max:",cf_max[period]), stringout("max_coeff.txt",s)/* save output to file as Maxima expressions */ );
[edit] More tutorials and code
- in Basic see Mandelbrot Dazibao
- in Java see Evgeny Demidov
- in C see Linas Vepstas
- in C++ see Wolf Jung page,
- in Factor program by Slava Pestov
- in Gnuplot see Tutorial by T.Kawano
- in Lisp (Maxima) see :
- in Octave see wikibooks or another verion by Christopher Wellons
- How to Plot the Mandelbrot Set By Hand
- comparison of various languages :
- rosettacode
- Fractal Benchmark in Ruby, Io, PHP, Python, Lua, Java, Perl, Applescript, TCL, ELisp, Javascript, OCaml, Ghostscript, and C by Erik Wrenholt
- in PDL, IDL, MATLAB, Octave, C and FORTRAN77 by Xavier Calbet
- ASCI graphic :
- using scripting languages by Warp
- using Lisp on Bill Clementson's blog
[edit] References
- ↑ Drawing the Mandelbrot set by the method of escape lines. M. Romera et al.
- ↑ http://www.metabit.org/~rfigura/figura-fractal/math.html boundary trace by Robert Figura
- ↑ http://www.geocities.com/CapeCanaveral/Launchpad/5113/fr27.htm| An open letter to Dr. Meech from Joyce Haslam in FRACTAL REPORT 27
- ↑ Heinz-Otto Peitgen (Editor, Contributor), Dietmar Saupe (Editor, Contributor), Yuval Fisher (Contributor), Michael McGuire (Contributor), Richard F. Voss (Contributor), Michael F. Barnsley (Contributor), Robert L. Devaney (Contributor), Benoit B. Mandelbrot (Contributor) : The Science of Fractal Images. Springer; 1 edition (July 19, 1988), page 199
- ↑ | The Mandelbrot Function by John J. G. Savard
- ↑ Wolf Jung explanation
- ↑ An algorithm to draw external rays of the Mandelbrot set by Tomoki Kawahira
- ↑ http://www.robertnz.net/cx.htm contour integration by Robert Davies
- ↑ [and exterior distance bounds for the Mandelbrot by Albert Lobo Cusido]
- ↑ Fractint : Misc. Options and algorithms
- ↑ Fractint doc by Noel Giffin
- ↑ Fractint doc by Noel Giffin
- ↑ A Series of spiral bays in the Mandelbrot set by Patrick Hahn
- ↑ Mark McClure "Bifurcation sets and critical curves" - Mathematica in Education and Research, Volume 11, issue 1 (2006).
- ↑ Burns A M : Plotting the Escape: An Animation of Parabolic Bifurcations in the Mandelbrot Set. Mathematics Magazine, Vol. 75, No. 2 (Apr., 2002), pp. 104-116
- ↑ Newsgroups: gmane.comp.mathematics.maxima.general Subject: system of equations Date: 2008-08-11 21:44:39 GMT
- ↑ Thayer Watkins : The Structure of the Mandelbrot Set
- ↑ Brown Method by Robert P. Munafo
- ↑ Fractal Stream help file
- ↑ [http://www.iec.csic.es/~gerardo/publica/Alvarez98a.pdf Determination of Mandelbrot Set's Hyperbolic Component Centres by G Alvarez, M Romera, G Pastor and F Montoya. Chaos, Solitons & Fractals 1998, Vol 9 No 12 pp 1997-2005]
- ↑ Myrberg, P J, Iteration der reelen polynome zweiten GradesIII, Ann. Acad. Sci. Fennicae, A. I no. 348, 1964.
- ↑ V Dolotin , A Morozow : On the shapes of elementary domains or why Mandelbrot set is made from almost ideal circles ?

using numerical methods




