Fractals/Iterations in the complex plane/MandelbrotSetExterior

From Wikibooks, open books for an open world
Jump to navigation Jump to search

Colouring of exterior of Mandelbrot set can be :

  • non-smooth :
    • Boolean/binary Escape Time Method ( bETM/M )
    • discrete = Level Set Method = LSM/M = iETM/M
  • Smooth :
    • radial measures
      • Real Escape Time Method( rETM/M )
      • Distance Estimation Method( DEM/M )
      • radius of complex potential = Complex Potential Method ( CPM/M )
  • angular measures
      • argument of complex potential
    • SAC = stripe average coloring
    • other

One can also draw curves :

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



Similar projects:


Escape time or dwell[edit]

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

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

Boolean escape time[edit]

This algorithm answers the question: “For which values of c will the Julia fractal, J(c), be line-like and for which dust-like?”[6]

Here complex plane consists of 2 sets : Mandelbrot set and its complement  :

ASCI graphic ( on screen)[edit]

ASCI graphic : Boolean escape time in text mode
// http://mrl.nyu.edu/~perlin/
main(k){float i,j,r,x,y=-16;while(puts(""),y++<15)for(x
=0;x++<84;putchar(" .:-;!/>)|&IH%*#"[k&15]))for(i=k=r=0;
j=r*r-i*i-2+x/25,i=2*r*i+y/10,j*j+i*i<11&&k++<111;r=j);}
-- 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 [7][8]

Graphic file ( PPM )[edit]

Here are various programs for creating pbm file [9]

C[edit]

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

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

===Integer escape time = LSM/M = dwell bands===
<gallery>
File:Mandelbrot seq.png|Number of details is proportional to maximal number of iterations
File:Animation of the growth of the Mandelbrot set as you iterate towards infinity.gif|Mandelbrot animation based on a static number of iterations per pixel. Here you can see why '''offset''' is sometimes used ( because - color gradient changes : for high MaxIteration disapears.
</gallery>
Here color is proportional to last iteration ( of final_n, final iteration).<ref>[http://plus.maths.org/content/computing-mandelbrot-set Computing the Mandelbrot set by Andrew Williams]</ref>

This is also called [[w:Level Set Method|Level Set Method]] ( LSM )

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

====C====

[[Image:Mandel_lsm_bw.jpg|thumb|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 :
<syntaxhighlight lang=C>
 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 C function whithout explicit complex numbers, only doubles:

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

here a short code with complex numbers:

// https://gitlab.com/adammajewski/mandelbrot_wiki_ACh/blob/master/betm.c
int iterate(double complex C , int iMax)
  {
   int i;
   double complex Z= 0.0; // initial value for iteration Z0
   
   for(i=0;i<iMax;i++)
    {
      Z=Z*Z+C; // https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
      if(cabs(Z)>EscapeRadius) break;
    }
   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);
}

A touch more optimised :

// optimised from cpp code by Geek3
inline int fractal(double cReal, double cImg, int max_iters)
// gives last iteration 
{
	double zReal = 0, zImg = 0, zReal2 = 0, zImg2 = 0;
	//iteration zero is always 0^2+0^2, it will never escape
	for (int it = 1; it < max_iters; it++)
	{	//because we have zReal^2 and zImg^2 pre-calculated
		//we can caclulate zImg first
		//then we don't need to calculate/store the "old" zReal
		zImg = (2 * zReal * zImg ) + cImg;
		zReal = zReal2 - zImg2 + cReal;

		// calculate next iteration: zReal^2 and zImg^2
		// they are used twice so calculate them once
		zReal2 = zReal * zReal;
		zImg2 = zImg * zImg;
		if (zReal2 + zImg2 > 4.0) return(it);
	}
	return(max_iters);
}
See also :

GLSL[edit]

Java[edit]

//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 Script[edit]

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.

Khan Academy[edit]

Lisp program[edit]

Whole Lisp program making ASCII 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 Gimp[edit]

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-Ray[edit]

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

Matemathica[edit]

Here is code by Paul Nylander

Level Curves of escape time Method = LCM/M[edit]

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 :

  • 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 , see explanation and source code. This method is very complicated for iterations > 5.

Decomposition of exterior of Mandelbrot set[edit]

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/M[edit]

binary decomposition: image with full code in C

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

also GLSL code from Fragmentarium :

#include "2D.frag"
#group Simple Mandelbrot

// maximal number of iterations
uniform int iMax; slider[1,100,1000] // increase iMax
// er2= er^2 wher er= escape radius = bailout
uniform float er2; slider[4.0,1000,10000] // increase er2

// compute color of pixel
 vec3 color(vec2 c) {
   vec2 z = vec2(0.0);  // initial value

 // iteration
  for (int i = 0; i < iMax; i++) {
    z = vec2(z.x*z.x-z.y*z.y,2*z.x*z.y) +  c; // z= z^2+c
    if (dot(z,z)> er2)   // escape test 
      // exterior
      if (z.x>0){ return vec3( 1.0);} // upper part of the target set 
      else return vec3(0.0); //lower part of the target set 
  }
  return vec3(0.0); //interior
}

Point c is plotting white or black if imaginary value of last iteration ( Zy) is positive or negative.[16]

nth-decomposition[edit]

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 Time[edit]

Other names of this method/algorithm are :

  • the fully-renormalized fractional iteration count ( by Linas Vepstas in 1997)[18]
  • smooth iteration count for generalized Mandelbrot sets ( by Inigo Quilez in 2016)[19]
  • continuous iteration count for the Mandelbrot set
  • Normalized Iteration Count Algorithm
  • Continuous coloring
  • smooth colour gradient
  • fractional iterations
  • fractional escape time

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

Other methods and speedups

Colouring formula in Ultrafractal :[20]

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

where :

  • log(log(bailout) can be precalculated

C[edit]

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[21]

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);
</source >

based on equation <ref>[http://linas.org/art-gallery/escape/escape.html Linas Vepstas. Renormalizing the mandelbrot escape.]</ref>
: <math>\nu(z) = n - \log_2 \log (z_n)\,</math>

====C++====
<syntaxhighlight lang=cpp>
// 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));

java[edit]

 /**
Smooth coloring algorithm
https://gitlab.com/shreyas.siravara/mandelbrot-with-smooth-coloring/blob/master/Mandelbrot.java
Mandelbrot with Smooth Coloring by Shreyas Siravara

*/
                double nsmooth = (iterations + 1 - Math.log(Math.log(Zn.getMagnitude())) / Math.log(ESCAPE_RADIUS));

                double smoothcolor = nsmooth / MAX_ITERATIONS;

                if (iterations < MAX_ITERATIONS) {
                    int rgb = Color.HSBtoRGB((float) (0.99f + 1.9 * smoothcolor), 0.9f, 0.9f);
                    g2d.setColor(new Color(rgb));
                } else {
                    g2d.setColor(Color.black.darker());
                }

Matemathica[edit]

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

Python[edit]

Python code using mpmath library[22]

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/M[edit]

Variants :

  • exterior DEM/M
  • interior DEM/M

Description

Complex potential[edit]

Description

See also[edit]

References[edit]

  1. Mathematics of Divergent Fractals by
  2. jussi harkonen : coloring-techniques
  3. wikipedia : Orbit trap
  4. Mandelbrot Orbit Trap Rendering! Programming How-To Video by DKM101
  5. Java program by Dieter Röß showing result of changing initial point of Mandelbrot iterations
  6. Julia fractals in PostScript by Kees van der Laan, EUROTEX 2012 & 6CM PROCEEDINGS 47
  7. Fractal Benchmark by Erik Wrenholt
  8. 12-minute Mandelbrot: fractals on a 50 year old IBM 1401 mainframe
  9. Computer Language Benchmarks Game
  10. example-code-from-presentation-ways-of-seeing-julia-sets by ed Burke
  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. linas.org : Renormalizing the Mandelbrot Escape
  19. I Quilez : mset_smooth
  20. fractalforums : What range/precision for fractional escape counts for Mandelbrot/Julia sets?
  21. Making Mandelbrot Set Movies by Tony Finch
  22. mpmath Python library