Last modified on 17 October 2014, at 20:16

Fractals/Iterations in the complex plane/Mandelbrot set

This book shows how to code different algorithms for drawing parameter plane [1] ( Mandelbrot set [2] ) for complex quadratic polynomial [3].

One can find different types of points / sets on parameter plane [4]


Exterior of Mandelbrot setEdit

Colouring of exterior of Mandelbrot set can be :

  • non-smooth :
    • boolean ( bETM/M )
    • discrete ( LSM/M = iETM/M)
  • Smooth  :
    • Real Escape Time ( rETM/M )
    • Distance estimation ( DEM/M )
    • Complex potential ( CPM/M )
    • "triangle inequality"[5]
    • Orbit trap [6]

One can also draw curves :

  • external rays
  • equipotential lines ( closed curves - quasi circles)



Escape timeEdit

Here for given point c on parameter plane one checks how critical point  z = 0.0 behaves on dynamical plane under forward iteration. If you change initial point ypu will get different result [7]

To draw given plane one needs to check/scan (all) its points. See here for more details ( optimisation) Read definitions first.

Boolean escape timeEdit

Here complex plane consists of 2 sets : Mandelbrot set M\, and its complement M^c\, :

\mathbb{C}= M \cup M^c

ASCI graphic ( on screen)Edit

ASCI graphic : Boolean escape time in text mode


-- Haskell code by Ochronus
-- http://blog.mostof.it/mandelbrot-set-in-ruby-and-haskell/
import Data.Complex
mandelbrot a = iterate (\z -> z^2 + a) a !! 500
main = mapM_ putStrLn [[if magnitude (mandelbrot (x :+ y)) < 2 then '*' else ' '
                           | x <- [-2, -1.9685 .. 0.5]]
                           | y <- [1, 0.95 .. -1]]


; common lisp
(loop for y from -1.5 to 1.5 by 0.05 do
      (loop for x from -2.5 to 0.5 by 0.025 do
		(let* ((c (complex x y)) ; parameter
                   	(z (complex 0 0))
                   	(iMax 20) ; maximal number of iterations
			(i 0)) ; iteration number
 
		(loop  	while (< i iMax ) do 
			(setq z (+ (* z z) c)) ; iteration
			(incf i)
			(when (> (abs z) 2) (return i)))
	   ; color of pixel
           (if (= i iMax) (princ (code-char 42)) ; inside M
                          (princ (code-char 32))))) ; outside M
      (format t "~%")) ; new line

Comparison programs in various languages [8]

Graphic file ( PPM )Edit

Here are various programs for creating pbm file [9]

CEdit

This is complete code of C one file 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 z_0 = z_{cr} = 0 \,

It uses some speed_improvement. Instead of checking :

sqrt(Zx2+Zy2)<ER

it checks :

(Zx2+Zy2)<ER2 // ER2 = ER*ER

It gives the same result but is faster.



 /* 
 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>
 #include <math.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 bytes) */
                        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;
 }

Integer escape time = LSM/MEdit

Here color is proportional to last iteration.[10]

This is also called Level Set Method ( LSM )

L_n= \{ c : z_n \in T ~~\mbox{and} ~~ z_k \notin T ~~\mbox{where}~~ k<n \}\,


CEdit

LSM/M image with full code in C

Difference 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 */    
             };


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

C++Edit

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 than pseudocode), but it can be coded in other way which can be executed faster .

Here is faster code :

// based on cpp code by Geek3
inline int fractal(double cx, double cy, int max_iters)
// gives last iteration 
{
	double zx = 0, zy = 0;	
	if (zx * zx + zy * zy > 4) return(0); // it=0
	for (int it = 1; it < max_iters; it++)
	{	double zx_old = zx;
		zx = zx * zx - zy * zy;
		zy = 2 * zx_old * zy;
		zx += cx;
		zy += cy;
		if (zx * zx + zy * zy > 4.0) return(it);
	}
	return(max_iters);
}

JavaEdit

//Java code by Josef Jelinek 
// http://java.rubikscube.info/
 
 int mandel(double cx, double cy) {
    double zx = 0.0, zy = 0.0;
    double zx2 = 0.0, zy2 = 0.0;
    int iter = 0;
    while (iter < iterMax && zx2 + zy2 < 4.0) {
      zy = 2.0 * zx * zy + cy;
      zx = zx2 - zy2 + cx;
      zx2 = zx * zx;
      zy2 = zy * zy;
      iter++;
    }
    return iter;
  }

Java ScriptEdit

Here is JavaScript function which does not give last iteration but LastIteration modulo maxCol. It makes colour cycling ( if maxCol < maxIt ).

function iterate(Cr,Ci) {
// JavaScript function by Evgeny Demidov
// http://www.ibiblio.org/e-notes/html5/fractals/mandelbrot.htm
  var I=0, R=0,  I2=0, R2=0,   n=0;
  if (R2+I2 > max) return 0;
  do  {  I=(R+R)*I+Ci;  R=R2-I2+Cr;  R2=R*R;  I2=I*I;  n++;
  } while ((R2+I2 < max) && (n < maxIt) );
  if (n == maxIt) return maxCol;  else return n % maxCol;
}

Above functions do not use explicit definition of complex number.

Lisp programEdit

Whole Lisp program making ASCI graphic based on code by Frank Buss [11] [12]


; common lisp
(loop for y from -1.5 to 1.5 by 0.1 do
      (loop for x from -2.5 to 0.5 by 0.04 do
            (let* ((i 0)
                   (z (complex x y))
                   (c z))
              (loop while (< (abs
                              (setq z (+ (* z z) c)))
                             2)
                    while (< (incf i) 32))
              (princ (code-char (+ i 32))))) ; ASCII chars <= 32 contains non-printing characters
      (format t "~%"))


MathMap plugin for GimpEdit

filter mandelbrot (gradient coloration)
    c=ri:(xy/xy:[X,X]*1.5-xy:[0.5,0]);
    z=ri:[0,0]; # initial value z0 = 0 
    # iteration of z
    iter=0;
    while abs(z)<2 && iter<31
    do
        z=z*z+c;  # z(n+1) = fc(zn)
        iter=iter+1
    end;
    coloration(iter/32) # color of pixel
end

Pov-RayEdit

Pov-Ray has a built-in function mandel[13]

MatemathicaEdit

Here is code by Paul Nylander

Level Curves of escape time Method = LCM/MEdit

edge detection of Level sets
Lemniscates of Mandelbrot set

Lemniscates are boundaries of Level Sets of escape time ( LSM/M ). They can be drawn using 2 methods:

  • edge detection of Level sets.
    • Algorithm described in paper by M. Romera et al.[14] This method is fast and allows looking for high iterations.
  • boundary trace[15]
  • drawing curves L_n(T)=\{c: abs(z_n)=ER  \}\,, see explanation and source code. This method is very complicated for iterations > 5.

Decomposition of exterior of Mandelbrot setEdit

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.

Binary decomposition of LSM/MEdit

binary decomposition: image with full code in C

Here target set T\, on dynamic plane is divided into 2 parts (binary decomposition = 2-decomposition ):

  • upper half ( white) T_u = \{ z : |z| > ER ~~\mbox{and}~~ Im (z) > 0\}\,
  • lower half (black) T_l = \{ z : |z| > ER ~~\mbox{and}~~ Im (z) \le 0 \}\,

Division of target set induces decomposition of level sets L_n\, into 2^{n+1}\, parts:

  • L_{n,u}  =\{ c:  |z_n| > ER ~~\mbox{and}~~ Im (z_n) > 0 \}\, which is colored white,
  • L_{n,l} = \{ c : |z_n| > ER ~~\mbox{and}~~ Im (z_n) \le 0 \}\, which is colored black.


External rays of angles (measured in turns):

angle = (k / 2^n ) ~~\mbox{mod }~1\,

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.[16]

nth-decompositionEdit

This method is extension of binary decomposition.

The target set T = { z : |zn| > R } with a very large escape radius ( for example R = 12 ) is divided into more then 2 parts ( for example 8).[17]

Real Escape TimeEdit

Other names of this method/algorithm are :

  • Normalized Iteration Count Algorithm
  • Continuous coloring
  • smooth colour gradient
  • fractional iterations

Here color of exterior of Mandelbrot set is proportional not to Last Iteration ( which is integer number) but to real number :

\nu(z) = \lim_{i \to \infty} (i - \log_2 \log_2 |z_i|)\,

Other methods and speedups

Colouring formula in Ultrafractal :[18]

smooth iter = iter + 1 + ( log(log(bailout)-log(log(cabs(z))) )/log(2)

where :

  • log(log(bailout) can be precalculated

CEdit

To use log2 function add :

#include <math.h>

at the beginning of program.


if (Iteration==IterationMax)
 { /*  interior of Mandelbrot set = black */
  color[0]=0;
  color[1]=0;
  color[2]=0;                           
 }
 /* exterior of Mandelbrot set  */
 else GiveRainbowColor((double)(Iteration- log2(log2(sqrt(Zx2+Zy2))))/IterationMax,color);

where :

  • Zx2 = Zx*Zx
  • Zy2 = Zy*Zy


Here is another version by Tony Finch[19]

while (n++ < max &&
 x2+y2 < inf) {
y = 2*x*y + b;
x = x2-y2 + a;
y2 = y*y;
x2 = x*x;
}
nu = n - log(log(x2+y2)/2)/ log(2);

based on equation [20]

\nu(z) = n - \log_2 \log (z_n)\,

C++Edit

// based on cpp code by Geek3 from http://en.wikibooks.org/wiki/File:Mandelbrot_set_rainbow_colors.png
sqrxy = x * x + y * y;
double m = LastIteration + 1.5 - log2(log2(sqrxy));

MatemathicaEdit

Here is code by Paul Nylander. It uses different formula :

cet = n + log_2ln(R) - log_2ln|z|

PythonEdit

Python code using mpmath library[21]

def mandelbrot(z):
    c = z
    for i in xrange(ITERATIONS):
        zprev = z
        z = z*z + c
        if abs(z) > ESCAPE_RADIUS:
            return ctx.exp(1j*(i + 1 - ctx.log(ctx.log(abs(z)))/ctx.log(2)))
    return 0

Distance estimation DEM/MEdit

Variants :

  • exterior DEM/M
  • interior DEM/M

Description

Complex potentialEdit

Uniformization of complement of Mandelbrot set

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 its absolute value. So on one image one can use more than one variable to color image.[22]

Real potential = CPM/MEdit

In Fractint :

potential =  log(modulus)/2^iterations

One can use real potential to:

  • smooth (continuous) coloring[23]
  • 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;

External angle and external ( parameter) rayEdit

First find angle of last iteration. It is easy to compute and shows some external rays as a borders of level sets.

Then one go futher.


MethodsEdit

TestsEdit

The Wolf Jung testEdit
Part of parameter plane with external rays 1/7, 321685687669320/2251799813685247 and 321685687669322/2251799813685247 landing on the Mandelbrot set

The external parameter rays for angles (in turns)

  • 321685687669320/2251799813685247 (period 51, lands on c1 = -0.088891642419446 +0.650955631292636i )
  • 321685687669322/2251799813685247 ( period 51 lands on c2 = -0.090588078906990 +0.655983860334813i )
  • 1/7 ( period 3, lands on c3 = -0.125000000000000 +0.649519052838329i )


Angles differ by about 10^{-15}, but the landing points of the corresponding parameter rays are about 0.035 apart. [27] It can be computed with Maxima CAS :

(%i1) c1: -0.088891642419446  +0.650955631292636*%i;
(%o1) 0.650955631292636*%i−0.088891642419446
(%i2) c2:-0.090588078906990  +0.655983860334813*%i;
(%o2) 0.655983860334813*%i−0.09058807890699
(%i3) abs(c2-c1);
(%o3) .005306692383854863
(%i4) c3: -0.125000000000000  +0.649519052838329*%i$
(%i5) abs(c3-c1);
(%o5) .03613692356607755
(%i6) a3:1/7$
(%i7) float(abs(a3-a1));
(%o7) 4.440892098500628*10^−16

Informations from W Jung program :

The angle  1/7  or  p001 has  preperiod = 0  and  period = 3.
The conjugate angle is  2/7  or  p010 .
The kneading sequence is  AA*  and the internal address is  1-3 .
The corresponding parameter rays are landing at the root of a satellite component of period 3.
It is bifurcating from period 1.
The angle  321685687669320/2251799813685247  or  p001001001001001001001001001001001001001001001001000 has  preperiod = 0  and  period = 51.
The conjugate angle is  321685687669319/2251799813685247  or  p001001001001001001001001001001001001001001001000111 .
The kneading sequence is  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABB*  and the internal address is  1-49-50-51 .
The corresponding parameter rays are landing at the root of a primitive component of period 51.
The angle  321685687669322/2251799813685247  or  p001001001001001001001001001001001001001001001001010 has  preperiod = 0  and  period = 51.
The conjugate angle is  321685687669329/2251799813685247  or  p001001001001001001001001001001001001001001001010001 .
The kneading sequence is  AABAABAABAABAABAABAABAABAABAABAABAABAABAABAABAABAA*  and the internal address is  1-3-51 .
The corresponding parameter rays are landing at the root of a satellite component of period 51.
It is bifurcating from period 3.


The test by G. Pastor and Miguel RomeraEdit

The external parameter rays for angles (in turns)

  • 6871947673/34359738367 ( period 35 )
  • 9162596898/34359738367 ( period 35 )

the central babies Mandelbrot sets of the cauliflowers located at -0.153756141 + 1.030383223i

(not that 34359738367 = 2^35 - 1)

test by M. Romera,1 G. Pastor, A. B. Orue,1 A. Martin, M.-F. Danca,and F. MontoyaEdit
Part of parameter plane with external 5 rays landing on the Mandelbrot set.png

G Pastor gave an example of external rays for which the resolution of the IEEE 754 is not sufficient [28]:

  • \theta_{267}^- = 0.((001)^{88}010)_2 = \frac{33877456965431938318210482471113262183356704085033125021829876006886584214655562}{237142198758023568227473377297792835283496928595231875152809132048206089502588927}
  • \theta_{267}^+ = 0.((001)^{87}010001)_2 = \frac{33877456965431938318210482471113262183356704085033125021829876006886584214655569}{237142198758023568227473377297792835283496928595231875152809132048206089502588927}
  • \theta_{3}^- = 0.(001)_2 = \frac{1}{7} =   0.(142857)_{10} ( period 3, lands on root point of period 3 component c3 = -0.125000000000000 +0.649519052838329i )
  • \theta_{268}^- = 0.((001)^{88}0001)_2 = \frac{67754913930863876636420964942226524366713408170066250043659752013773168429311121}{474284397516047136454946754595585670566993857190463750305618264096412179005177855}
  • \theta_{268}^+ = 0.((001)^{88}0010)_2 = \frac{67754913930863876636420964942226524366713408170066250043659752013773168429311122}{474284397516047136454946754595585670566993857190463750305618264096412179005177855}


One can analyze these angles using program by Claude Heiland-Allen :

./bin/mandelbrot_describe_external_angle ".(001)"
binary: .(001)
decimal: 1/7
preperiod: 0
period: 3
 
 
 
./bin/mandelbrot_describe_external_angle 
".(001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001010)"
binary: 
.(001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001010)
decimal: 
33877456965431938318210482471113262183356704085033125021829876006886584214655562/237142198758023568227473377297792835283496928595231875152809132048206089502588927
preperiod: 0
period: 267
 
 
./bin/mandelbrot_describe_external_angle 
".(001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001010001)"
binary: 
.(001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001001010001)
decimal: 
33877456965431938318210482471113262183356704085033125021829876006886584214655569/237142198758023568227473377297792835283496928595231875152809132048206089502588927
preperiod: 0
period: 267
 
 
 
./bin/mandelbrot_describe_external_angle 
".(0010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010001)"
binary: 
.(0010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010001)
decimal: 
67754913930863876636420964942226524366713408170066250043659752013773168429311121/474284397516047136454946754595585670566993857190463750305618264096412179005177855
preperiod: 0
period: 268
 
 
 
./bin/mandelbrot_describe_external_angle 
".(0010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010)"
binary: 
.(0010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010010)
decimal: 
67754913930863876636420964942226524366713408170066250043659752013773168429311122/474284397516047136454946754595585670566993857190463750305618264096412179005177855
preperiod: 0
period: 268


Landing points of above rays are roots with angled internal addresses ( description by Claude Heiland-Allen) :

  • the upper one will be 1 -> 1/3 -> 3 -> 1/(period/3) -> period because it's the nearest bulb to the lower root cusp of 1/3 bulb and child bulbs of 1/3 bulb have periods 3 * denominator(internal angle) ie, 1 -> 1/3 -> 3 -> 1/89 -> 267
  • the lower one will be 1 -> floor(period/3)/period -> period because it's the nearest bulb below the 1/3 cusp ie, 1 -> 89/268 -> 268
  • the middle ray .(001) lands at the root of 1 -> 1/3 -> 3, from the cusp on the lower side (which is on the right in a standard unrotated view)

Boundary of Mandelbrot setEdit

description

Interior of Mandelbrot set - hyperbolic componentsEdit

The Lyapunov ExponentEdit

Lyapunov exponents of mini the Mandelbrot set
Lyapunov exponent of real quadratic map

Math equation :[29]

\lambda_f(z_0) = \lim_{n\rightarrow\infty}  \frac{1}{n} \sum_{i=0}^{n-1} \left ( \ln\left|f'(z_i)\right | \right )

where :

f'(x) = \frac{d}{dz}f_c(z) = 2z

means first derivative of f with respect to z


See also :

  • image and description by janthor [30]
  • image by Anders Sandberg [31]

Interior distance estimationEdit

Interior distance estimation

Description of method

absolute value of the orbitEdit

# Hypercomputing the Mandelbrot Set? by Petrus H. Potgieter February 1, 2008
n=1000; # For an nxn grid
m=50; # Number of iterations
c=meshgrid(linspace(-2,2,n))\ # Set up grid
+i*meshgrid(linspace(2,-2,n));
x=zeros(n,n); # Initial value on grid
for i=1:m
x=x.^2+c; # Iterate the mapping
endfor
imagesc(min(abs(x),2.1)) # Plot monochrome, absolute
# value of 2.1 is escape


internal level setsEdit

Color of point  :

  • is proportional to the value of z is at final iteration.
  • shows internal level sets of periodic attractors.


bof60Edit

Image of bof60 in on page 60 in the book "the Beauty Of Fractals".Description of the method described on page 63 of bof. It is used only for interior points of the Mandelbrot set.

Color of point is proportional to :

  • the smallest distance of its orbit from origin[32][33]
  • the smallest value z gets during iteration [34]
  • illuminating the closest approach the iterates of the origin (critical point) make to the origin inside the set


Level sets of distance are sets of points with the same distance[35]

if (Iteration==IterationMax)
 /* interior of Mandelbrot set = color is proportional to modulus of last iteration */
 else { /* exterior of Mandelbrot set = black */
  color[0]=0;
  color[1]=0;
  color[2]=0;                           
 }
  • fragment of code : fractint.cfrm from Gnofract4d [36]
bof60 {
 init:
       float mag_of_closest_point = 1e100
 loop:
       float zmag = |z|
       if zmag < mag_of_closest_point
               mag_of_closest_point = zmag
       endif
 final:
       #index = sqrt(mag_of_closest_point) * 75.0/256.0
}

bof61Edit

3D view of bof61
2D view of BOF61

This is the method described in the book "The Beauty of Fractals" on page 63, but the image in on page 61.

Color of point is proportional to :

  • the time it takes z to reach its smallest value
  • iterate of the critical point makes the closest approach
  • 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)[37] [38].


Fragment of code : fractint.cfrm from Gnofract4d [39]

bof61 {
 init:
       int current_index = -1 ; -1 to match Fractint's notion of iter count
       int index_of_closest_point = -1
       float mag_of_closest_point = 1e100
 loop:
       current_index = current_index + 1
       float zmag = |z|
       if zmag < mag_of_closest_point
               index_of_closest_point = current_index
               mag_of_closest_point = zmag
       endif
 final:
       #index = index_of_closest_point /256.0
}

Cpp function

//       function is based on function mclosetime
//       from mbrot.cpp 
//       from program mandel by Wolf Jung 
//      http:www.iram.rwth-aachen.de/~jung/indexp.html  
////8 = iterate =  bof61
// bailout2=4 
int mclosetime(std::complex<double> C , int iter_max,  int bailout2)
{ int j, cln = 0;
  double x = C.real(), y = C.imag(), u, cld = 4;
  for (j = 0; j <= iter_max; j++)
  { u = x*x + y*y; 
    if (u > bailout2) return j;
    if (u < cld) {cld = u;cln = j; }
 
    u = x*x - y*y + C.real();
    y = 2*x*y + C.imag();
    x = u;
  }
  return iter_max + cln % 15; //iterate, bof61
 
}

It can be used :

// compute escape time  
              int last_iter= mclosetime(C,iter_max,bailout2);
              //  drawing code */
              if (last_iter>=iter_max) { putpixel(ix,iy,last_iter - iter_max);} // interior
                                else putpixel(ix,iy,WHITE); // exterior

Period of hyperbolic componentsEdit

period of hyperbolic components

Period of hyperbolic component of Mandelbrot set is a period of limit set of critical orbit.


Algorithms for computing period:

  • direct period detection from iterations of critical point z = 0.0 on dynamical plane
  • "quick and dirty" algorithm : check if abs(z_n ) < eps   then colour c-point with colour n. Here n is a period of attracting orbit and eps is a radius of circle around attracting point = precision of numerical computations
  • "methods based on interval arithmetic when implemented properly are capable of finding all period-n cycles for considerable large n." (ZBIGNIEW GALIAS )[40]
  • Floyd's cycle-finding algorithm [41]
  • the spider algorithm

Multiplier mapEdit

definition

Internal angleEdit

Method by Renato Fonseca : [42] "a point c in the set is given a hue equal to argument

arg(z_{n_{max}}) = arctan\frac{Im(z_{n_{max}})}{Re(z_{n_{max}})}

(scaled appropriatly so that we end up with a number in the range 0 - 255). The number z_nmax is the last one calculated in the z's sequence. "

Internal raysEdit

Internal and external rays

When radius\, varies and angle\, is constant then c\, goes along internal ray. It is used as a path inside Mandelbrot set

/* find c in component of Mandelbrot set 
 uses complex type so #include <complex.h> and -lm 
 uses code by Wolf Jung from program Mandel
 see function mndlbrot::bifurcate from mandelbrot.cpp
 http://www.mndynamics.com/indexp.html
 
  */
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int period)
{
  //0 <= InternalRay<= 1
  //0 <= InternalAngleInTurns <=1
  double t = InternalAngleInTurns *2*M_PI; // from turns to radians
  double R2 = InternalRadius * InternalRadius;
  double Cx, Cy; /* C = Cx+Cy*i */
  switch ( period ) {
    case 1: // main cardioid
      Cx = (cos(t)*InternalRadius)/2-(cos(2*t)*R2)/4; 
      Cy = (sin(t)*InternalRadius)/2-(sin(2*t)*R2)/4; 
      break;
   case 2: // only one component 
      Cx = InternalRadius * 0.25*cos(t) - 1.0;
      Cy = InternalRadius * 0.25*sin(t); 
      break;
  // for each period  there are 2^(period-1) roots. 
  default: // safe values
      Cx = 0.0;
      Cy = 0.0; 
    break; }
 
  return Cx+ Cy*I;
}
 
// draws points to memmory array data
int DrawInternalRay(double InternalAngleInTurns , unsigned int period, int iMax, unsigned char data[])
{
 
   complex double c;
   double InternalRadius;
   double RadiusStep; // between radius of points 
   int i; // number of point to draw
 
  RadiusStep = 1.0/iMax;
 
  for(i=0;i<=iMax;++i){ 
   InternalRadius = i * RadiusStep;
   c = GiveC(InternalAngleInTurns, InternalRadius, period);
   DrawPoint(c,data);
  }
 
return 0;
}

Example : internal ray of angle =1/6 of main cardioid.

Internal angle :

angle = 1/6 \,

radius of ray :

 0 \le radius \le 1  \,

Point of internal radius of unit circle :

 w = radius * e^{i * angle}\,

Map point w to parameter plane :

c = \frac{w}{2} - \frac{w^2}{4} \,

For epsilon = 0 \, this is equation for main cardioid.

Internal curveEdit

When radius\, is constant varies and angle\, varies then c\, goes along internal curve.

/* find c in component of Mandelbrot set 
 uses complex type so #include <complex.h> and -lm 
 uses code by Wolf Jung from program Mandel
 see function mndlbrot::bifurcate from mandelbrot.cpp
 http://www.mndynamics.com/indexp.html
 
  */
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int period)
{
  //0 <= InternalRay<= 1
  //0 <= InternalAngleInTurns <=1
  double t = InternalAngleInTurns *2*M_PI; // from turns to radians
  double R2 = InternalRadius * InternalRadius;
  double Cx, Cy; /* C = Cx+Cy*i */
  switch ( period ) {
    case 1: // main cardioid
      Cx = (cos(t)*InternalRadius)/2-(cos(2*t)*R2)/4; 
      Cy = (sin(t)*InternalRadius)/2-(sin(2*t)*R2)/4; 
      break;
   case 2: // only one component 
      Cx = InternalRadius * 0.25*cos(t) - 1.0;
      Cy = InternalRadius * 0.25*sin(t); 
      break;
  // for each period  there are 2^(period-1) roots. 
  default: // safe values
      Cx = 0.0;
      Cy = 0.0; 
    break; }
 
  return Cx+ Cy*I;
}
 
 
// draws points to memmory array data
int DrawInternalCurve(double InternalRadius , unsigned int period,  int iMax, unsigned char data[])
{
 
   complex double c;
   double InternalAngle; // in turns = from 0.0 to 1.0
   double AngleStep;
   int i;
   // int iMax =100;
 
  AngleStep = 1.0/iMax;
 
  for(i=0;i<=iMax;++i){ 
   InternalAngle = i * AngleStep;
   c = GiveC(InternalAngle, InternalRadius, period);
   DrawPoint(c,data);
  }
 
return 0;
}

Boundaries of hyperbolic componentsEdit

description

Centers of componentsEdit

Methods of finding centers :


  • all centers for given period p\,
    • solving F_p(0,c) = 0 \, or   G_p(0,c) = 0\,
      • showing results in graphical form
      • using numerical methods[43]
        • "quick and dirty" algorithm : check if abs(z_n ) < eps then colour c-point with colour n. Here n is a period of attracting orbit and eps is a radius of circle around attracting point = precision of numerical comutations
        • finding all roots of polynomial
        • Myrberg's method [44][45]
      • interval arithmetics and a divide and conquere strategy[46]
  • one center c_p\, for given period p\, near given point c\,


System of equations defining centersEdit

Centers of hyperbolic components for period 10
Centers of hyperbolic components for period 1-10 computed using irreducible polynomials

Center of hyperbolic component is a point c\, of parameter plain, for which periodic z-cycle is superattracting. It gives a system of 2 equations defining centers of period p\, hyperbolic components :

  • first defines periodic point z_0\,,
  • second makes z_0\, point superattracting.

\begin{cases}
    F_p(z_0,c) = z_0 \\
     \frac{df_c^{(p)}(z_0)}{dz}=0 
\end{cases}


see definitions for details.

Solving system of equationsEdit

Second equation can be solved when critical point z_{cr} = 0 \, is in this cycle :

z_0 = z_{cr} = 0 \,

To solve system put z_0\, into first equation.

Equation defining centersEdit

One gets equation :

F_p(0,c) = 0 \,

Centers of components are roots of above equation.

Because z = 0 \, one can remove z from these equations :

  • for period 1 : z^2+c=z and z=0 so c=0
  • for period 2 : (z^2+c)^2 +c =z and z=0 so c^2 + c =0
  • for period 3 : ((z^2+c)^2 +c)^2 +c = z and z=0 so (c^2 + c)^2 +c =0

Here is Maxima functions for computing above functions :

P[n]:=if n=1 then c else P[n-1]^2+c;
Reduced equation defining centersEdit

Set of roots of above equation contains centers for period p and its divisors. It is because :

F_p(0,c) = \prod_{m|p} G_m(0,c) \,

where :

For example :


F_1 = c = G_1 \,

F_2 = c*(c+1) = G_1* G_2 \,

F_3 = c*(c^3+2*c^2+c+1) = G_1* G_3 \,

F_4 = c*(c+1)*(c^6+3*c^5+3*c^4+3*c^3+2*c^2+1) = G_1* G_2* G_4 \,

So one can find irreducible polynomials using :


G_p(0,c) = \frac{F_p(0,c) }{ \prod_{m|p,m<p} G_m(0,c)}  \,


Here is Maxima function for computing G \, :


GiveG[p]:=
block(
[f:divisors(p),
t:1], /* t is temporary variable = product of Gn for (divisors of p) other than 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 F\, and G\, for periods 1-10 and precision needed to compute roots of these functions.


period p \, \deg(F_p(0,c))  \, \deg(G_p(0,c))  \, fpprec \,
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 p \, log10(gc_p)\, gc_p \,
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 than size of binary form of greatest coefficient log_2(gc)\, :

fpprec > log_2(gc)\,
period p \, \deg(G_p(0,c))  \, \log_2(gc_p)  \, fpprec \,
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 gc\,

period_Max:11;
/* ----------------- definitions -------------------------------*/
/* basic function  = 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 */
);

Graphical methods for finding centersEdit

All these methods shows centers for period n and its divisors.

Animation showing colors proportional to absolute value of iterations 1-20
Absolute value of first iteration
Color shows in which quadrant 5-th iteration lands
color is proportional to magnitude of znEdit
  • color is proportional to magnitude of zn[49]
  • Parameter plane is scanned for points c for which orbit of critical point vanishes [50]
  • youtube video : Mandelbrot Oscillations [51]
color shows in which quadrant zn landsEdit
ninth-decomposition. Image and src code

This is radial nth-decomposition of exterior of Mandelbrot set ( compare it with n-th decomposition of LSM/M )

4 colors are used because there are 4 quadrants :

  • re(z_n) > 0 and im(z_n) > 0 ( 1-st quadrant )
  • re(z_n) < 0 and im(z_n) > 0 ( 2-nd quadrant )
  • re(z_n) < 0 and im(z_n) < 0 ( 3-rd quadrant )
  • re(z_n) > 0 and im(z_n) < 0 (4-th quadrant ).


"... when the four colors are meeting somewhere, you have a zero of q_n(c), i.e., a center of period dividing n. In addition, the light or dark color shows if c is in M." ( Wolf Jung )


Here is fragment of c code for computing 8-bit color for point c = cx + cy * i  :

/* initial value of orbit = critical point Z= 0 */
                       Zx=0.0;
                       Zy=0.0;
                       Zx2=Zx*Zx;
                       Zy2=Zy*Zy;
                       /* the same number of iterations for all points !!! */
                       for (Iteration=0;Iteration<IterationMax; Iteration++)
                       {
                           Zy=2*Zx*Zy + Cy;
                           Zx=Zx2-Zy2 +Cx;
                           Zx2=Zx*Zx;
                           Zy2=Zy*Zy;
                       };
                       /* compute  pixel color (8 bit = 1 byte) */
                       if ((Zx>0) && (Zy>0)) color[0]=0;   /* 1-st quadrant */
                       if ((Zx<0) && (Zy>0)) color[0]=10; /* 2-nd quadrant */
                       if ((Zx<0) && (Zy<0)) color[0]=20; /* 3-rd quadrant */
                       if ((Zx>0) && (Zy<0)) color[0]=30; /* 4-th quadrant */

One can also find cpp code in function quadrantqn from class mndlbrot in mndlbrot.cpp file from source code of program Mandel by Wolf Jung[52]

The differences from standard loop are :

  • there is no bailout test (for example for max value of abs(zn) ). It means that every point has the same number of iterations. For large n overflow is possible. Also one can't use test Iteration==IterationMax
  • Iteration limit is relatively small ( start for example IterationMax = 8 )

See also code in Sage[53]

using Gaussian wave function method [54]Edit

VideoEdit

One can make videos using :

  • going along some paths on parameter plane ( for example internal and external rays )
  • zoom into parameter plane[55][56][57] using automatic determination of Iteration Max number[58]
  • changing coloring scheme ( for example color cycling - Fractint)
  • changing some parameters of algorithm, for example :
    • maximal iteration of escape time algorithm
    • bailout value [59]

Speed improvements - optimisationEdit

SymmetryEdit

The Mandelbrot set is symmetric with respect to the x-axis in the plane :

colour(x,y) = colour(x,-y)

its intersection with the x-axis ( real slice of Mandelbrot set ) is an interval :

<-2 ; 1/4>

It can be used to speed up computations ( up to 2-times )[60]

Bailout testEdit

Instead of checking :

 \sqrt{Z_x^2 + Z_y^2} < ER 

compute ER2 once and check :

 Z_x^2 + Z_y^2 < ER^2 

It gives the same result and is faster.

Period detectionEdit

"When rendering a Mandelbrot or Julia set, the most time-consuming parts of the image are the “black areas”. In these areas, the iterations never escape to infinity, so every pixel must be iterated to the maximum limit. Therefore, much time can be saved by using an algorithm to detect these areas in advance, so that they can be immediately coloured black, rather than rendering them in the normal way, pixel by pixel. FractalNet uses a recursive algorithm to split the image up into blocks, and tests each block to see whether it lies inside a “black area”. In this way, large areas of the image can be quickly coloured black, often saving a lot of rendering time. ... (some) blocks were detected as “black areas” and coloured black immediately, without having to be rendered. (Other) blocks were rendered in the normal way, pixel by pixel." Michael Hogg [61]


 // cpp code by Geek3 
// http://commons.wikimedia.org/wiki/File:Mandelbrot_set_rainbow_colors.png
bool outcircle(double center_x, double center_y, double r, double x, double y)
{ // checks if (x,y) is outside the circle around (center_x,center_y) with radius r
        x -= center_x;
        y -= center_y;
        if (x * x + y * y > r * r)
                return(true);
        return(false);
 
 // skip values we know they are inside
                        if ((outcircle(-0.11, 0.0, 0.63, x0, y0) || x0 > 0.1)
                                && outcircle(-1.0, 0.0, 0.25, x0, y0)
                                && outcircle(-0.125, 0.744, 0.092, x0, y0)
                                && outcircle(-1.308, 0.0, 0.058, x0, y0)
                                && outcircle(0.0, 0.25, 0.35, x0, y0))
                        {
                          // code for iteration
                         }

Cardioid and period-2 checkingEdit

One way to improve calculations is to find out beforehand whether the given point lies within the cardioid or in the period-2 bulb. Before passing the complex value through the escape time algorithm, first check if:

 (x+1)^2 + y^2 < \frac{1}{16}

to determine if the point lies within the period-2 bulb and

 q = \left(x - \frac{1}{4}\right)^2 + y^2
 q \left(q + \left(x - \frac{1}{4}\right)\right) < \frac{1}{4}y^2.

to determine if the point lies inside the main cardioid.

Periodicity checkingEdit

Most points inside the Mandelbrot set oscillate within a fixed orbit. There could be anything from ten to tens of thousands of points in between, but if an orbit ever reaches a point where it has been before then it will continually follow this path, will never tend towards infinity and hence is in the Mandelbrot set. This Mandelbrot processor includes periodicity checking (and p-2 bulb/cardioid checking) for a great speed up during deep zoom animations with a high maximum iteration value.


Perturbation theoryEdit

Very highly magnified images require more than the standard 64-128 or so bits of precision most hardware floating-point units provide, requiring renderers use slow "bignum" or "arbitrary precision"[62] math libraries to calculate. However, this can be sped up by the exploitation of perturbation theory[63]. Given

 z_{n+1} = z_n^2 + c

as the iteration, and a small epsilon, it is the case that

 (z_n + \epsilon)^2 + c = z_n^2 + 2z_n\epsilon + \epsilon^2 + c

or

 z_{n+1} + 2z_n\epsilon + \epsilon^2

so if one defines

 \epsilon_{n+1} = 2z_n\epsilon_n + \epsilon_n^2

one can calculate a single point (e.g. the center of an image) using normal, high-precision arithmetic (z), giving a reference orbit, and then compute many points around it in terms of various initial offsets epsilon-zero plus the above iteration for epsilon. For most iterations, epsilon does not need more than 16 significant figures, and consequently hardware floating-point may be used to get a mostly accurate image.[64] There will often be some areas where the orbits of points diverge enough from the reference orbit that extra precision is needed on those points, or else additional local high-precision-calculated reference orbits are needed. This rendering method, and particularly the automated detection of the need for additional reference orbits and automated optimal selection of same, is an area of ongoing, active research. Renderers implementing the technique are publicly available and offer speedups for highly magnified images in the multiple orders of magnitude range.[65]

More tutorials and codeEdit

ReferencesEdit

  1. parameter plane in wikipedia
  2. Mandelbrot set in wikipedia
  3. complex quadratic polynomial in wikipedia
  4. reenigne blog : mandelbrot-set-taxonomy
  5. Mathematics of Divergent Fractals by
  6. wikipedia : Orbit trap
  7. Java program by Dieter Röß showing result of changing initial point of Mandelbrot iterations
  8. Fractal Benchmark by Erik Wrenholt
  9. The Computer Language Benchmarks Game
  10. Computing the Mandelbrot set by Andrew Williams
  11. LIsp Program by Frank Buss
  12. Mandelbrot Set ASCII art at Bill Clementson's blog
  13. mandel function from 2.5.11.14 Fractal Patterns at Pov-Ray docs
  14. Drawing the Mandelbrot set by the method of escape lines. M. Romera et al.
  15. http://www.metabit.org/~rfigura/figura-fractal/math.html boundary trace by Robert Figura
  16. http://web.archive.org/20010415125044/www.geocities.com/CapeCanaveral/Launchpad/5113/fr27.htm%7C An open letter to Dr. Meech from Joyce Haslam in FRACTAL REPORT 27
  17. mandelbrot set n-th-decomposition
  18. fractalforums : What range/precision for fractional escape counts for Mandelbrot/Julia sets?
  19. Making Mandelbrot Set Movies by Tony Finch
  20. Linas Vepstas. Renormalizing the mandelbrot escape.
  21. mpmath Python library
  22. | The Mandelbrot Function by John J. G. Savard
  23. Smooth colouring is the key to the Mandelbrot set by Tony Finch
  24. An algorithm to draw external rays of the Mandelbrot set by Tomoki Kawahira
  25. Wolf Jung explanation
  26. [http://www.hindawi.com/journals/mpe/2013/105283/ A Method to Solve the Limitations in Drawing External Rays of the Mandelbrot Set M. Romera,1 G. Pastor, A. B. Orue,1 A. Martin, M.-F. Danca,and F. Montoya ]
  27. Wolf Jung's test for precision of drawing parameter rays
  28. A Method to Solve the Limitations in Drawing External Rays of the Mandelbrot Set M. Romera,1 G. Pastor, A. B. Orue,1 A. Martin, M.-F. Danca,and F. Montoya
  29. The logistic equation by Didier Gonze October 4, 2013
  30. Ljapunov Exponent and mandelbrot set by janthor
  31. Image by Anders Sandberg
  32. Fractint : Misc. Options and algorithms
  33. Java™ Number Cruncher: The Java Programmer's Guide to Numerical Computing By Ronald Mak
  34. Firefly Application Help by Terry W. Gintz
  35. Fractint doc by Noel Giffin
  36. gnofract4d
  37. Fractint doc by Noel Giffin
  38. A Series of spiral bays in the Mandelbrot set by Patrick Hahn
  39. gnofract4d
  40. Rigorous Investigations Of Periodic Orbits In An Electronic Circuit By Means Of Interval Methods by Zbigniew Galias
  41. Mandelbrot set drawing by Milan
  42. The Mandelbrot set by Renato Fonseca
  43. Piers Lawrence and Rob Corless : Mandelbrot Polynomials and Matrices. 1 ; 048 ; 575 roots of period 21
  44. [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]
  45. Myrberg, P J, Iteration der reelen polynome zweiten GradesIII, Ann. Acad. Sci. Fennicae, A. I no. 348, 1964.
  46. evaldraw script that can compute 1 million root by knighty
  47. irreducible divisors at wikipedia
  48. V Dolotin , A Morozow : On the shapes of elementary domains or why Mandelbrot set is made from almost ideal circles ?
  49. [http://www.cerebralmastication.com/2009/02/mandelbrot-set-in-r/ Mandelbrot Set in R by J.D. Long]
  50. Fractal Stream help file
  51. youtube video : Mandelbrot Oscillations
  52. Multiplatform cpp program Mandel by Wolf Jung
  53. Formation of Escape-Time Fractals By Christopher Olah
  54. [http://www.dhushara.com/DarkHeart/DarkHeart.htm#_3__The_Dark Exploding the Dark Heart of Chaos by Chris King March-Dec 2009]
  55. Really Deep Fractal Zoom Movie – Much Faster by Bruce Dawson
  56. Making Mandelbrot Set Movies by Tony Finch
  57. MLbrot by Daniel de Rauglaud
  58. Discussion : A way to determine the ideal number of maximum iterations for an arbitrary zoom level in a Mandelbrot fractal
  59. Gif image by jgabase : a wormhole effect on your fractals by changing the bailout dynamicaly
  60. How to use symetry of set
  61. FractalNet by Michael Hogg
  62. arbitrary precision at wikipedia
  63. perturbation theory in wikipedia
  64. "Superfractalthing - Arbitrary Precision Mandelbrot Set Rendering in Java". http://www.fractalforums.com/announcements-and-news/superfractalthing-arbitrary-precision-mandelbrot-set-rendering-in-java/. 
  65. "Kalles Fraktaler 2". http://www.chillheimer.de/kallesfraktaler/. 
  66. ASCII graphic

NotesEdit

  1. ^ It is conjectured that the polynomials (for the centers) are irreducible, but this has not been proved yet - Wolf Jung