Fractals/Iterations in the complex plane/demm

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

This algorithm has 2 versions :

  • exterior
  • interior

Compare it with version for dynamic plane and Julia set : DEM/J

Exterior DEM/M[edit]

Exterior DEM/M
simple boundary with DEM/M
Boundary with DEM/M and Sobel filter

Distance Estimation Method for Mandelbrots set ( DEM/M ) estimates distance from point c\, ( in the exterior of Mandelbrot set ) to nearest point in Mandelbrot set.

It can be used to create B&W images from BOF :[1]

  • map 41 on page 84
  • map 43 on page 85
  • an unnumbered plot on page 188

Math formule[edit]

Math formula : [2][3] [4]

distance(c_0,\sigma M) =\lim_{n \to \infty} 2 \frac{|z_n|}{|z'_n|}\ln|z_n|

where :

z_n = f_c^n(z=0.0)
z_n' = \frac{d}{dc} f_c^n(z_0).

is a first derivative of f_c^n(z_0) with respect to c which can be found by iteration starting with

z_0' = \frac{d}{dc} f_c^0(z_0) = 1

and then replacing at every consecutive step

z_{n+1}' = \frac{d}{dc} f_c^{n+1}(z_0) = 2\cdot{}f_c^n(z)\cdot\frac{d}{dc} f_c^n(z_0) + 1 = 2 \cdot z_n \cdot z_n' +1.


Pseudocode[edit]

  • The Beauty of Fractals
  • The Science of Fractal Images, page 198,
  • Distance Estimator by Robert P. Munafo[5]

Analysis of code from BOF[edit]

"The Beauty of Fractals gives an almost correct computer program for the distance estimation shown in the right image. A possible reason that that method did not gain ground is that the procedure in this program is seriously flawed: The calculation of z_k is performed (and completed) before the calculation of z'_k, and not after as it ought to be (z'_{k+1} uses z_k, not z_{k+1}). For the successive calculation of z'_k, we must know f'(z_k) (which in this case is 2z_k). In order to avoid the calculation of z_k (k = 0, 1, 2, ...) again, this sequence is saved in an array. Using this array, z'_{k+1} = 2z_kz'_k + 1 is calculated up to the last iteration number, and it is stated that overflow can occur. If overflow occurs the point is regarded as belonging to the boundary (the bail-out condition). If overflow does not occur, the calculation of the distance can be performed. Apart from it being untrue that overflow can occur, the method makes use of an unnecessary storing and repetition of the iteration, making it unnecessarily slower and less attractive. The following remark in the book is nor inviting either: "It turns out that the images depend very sensitively on the various choices" (bail-out radius, maximum iteration number, overflow, thickness of the boundary and blow-up factor). Is it this nonsense that has got people to lose all desire for using and generalizing the method? " Gertbuschmann

So :

  • in the loop, derivative should be calculated before the new z

Example code[edit]

Two algorithm in two loops[edit]


Here is Java function. It computes in one loop both : iteration z_n\, and derivative z'_n\,. If (on dynamic plane) critical point :

  • does not escapes to infinity ( bailouts), then it belongs to interior of Mandelbrot set and has colour maxColor,
  • escapes then it maybe in exterior or near boundary. Its colour is "proportional" to distance between the point c and the nearest point in the Mandelbrot set. It uses also colour cycling ( (int)R % maxColor ).
// Java function by Evgeny Demidov from http://www.ibiblio.org/e-notes/MSet/DEstim.htm
// based on code by Robert P. Munafo from http://www.mrob.com/pub/muency/distanceestimator.html
  public int iterate(double cr, double ci, double K, double f) {
    double Cr,Ci, I=0, R=0,  I2=I*I, R2=R*R, Dr=0,Di=0,D;   int n=0;
    if (f == 0){ Cr = cr; Ci = ci;}
    else{ Cr = ci; Ci = cr;}
    do  {
      D = 2*(R*Dr - I*Di) +1;  Di = 2*(R*Di + I*Dr);  Dr = D;
      I=(R+R)*I+Ci;  R=R2-I2+Cr;  R2=R*R;  I2=I*I;  n++;
    } while ((R2+I2 < 100.) && (n < MaxIt) );
    if (n == MaxIt) return maxColor; // interior
    else{ // boundary and exterior 
     R = -K*Math.log(Math.log(R2+I2)*Math.sqrt((R2+I2)/(Dr*Dr+Di*Di))); // compute distance
     if (R < 0) R=0;
     return (int)R % maxColor; }; 
 }


Here is cpp function. It gives integer index of color as an output.

/*
 this function is  from program mandel ver 5.10 by Wolf Jung
 see file mndlbrot.cpp   
 http://www.mndynamics.com/indexp.html
 
It is checked first (in pixcolor)
that the point escapes before calling this function.  
So we do not compute the derivative and the logarithm 
when it is not escaping.  
This would not be a good idea if most points escape anyway.
 
*/
int mndlbrot::dist(double a, double b, double x, double y)
{  uint j; 
   double xp = 1, yp = 0; // zp = xp+yp*i =  1 ; derivative 
   double nz, nzp; 
 
  // iteration 
  for (j = 1; j <= maxiter; j++)
   {  // zp
      nz = 2*(x*xp - y*yp); 
      yp = 2*(x*yp + y*xp); 
      xp = nz; //zp = 2*z*zp; on the dynamic plane 
      // if sign is positive it is parameter plane, if negative it is dynamic plane.
      if (sign > 0) xp++; //zp = 2*z*zp + 1 ; on the parameter plane
      //z = z*z + c; 
      nz = x*x - y*y + a; 
      y = 2*x*y + b; 
      x = nz; 
      //
      nz = x*x + y*y; 
      nzp = xp*xp + yp*yp;
      // 2 conditions for stopping the iterations
      if (nzp > 1e40 || nz > bailout) break;
   }
 
   // 5 types of points but 3 colors  
   /*  not escaping, rare 
       Is it not simply interior ???
       This should not be necessary.  I do not know why I included it,
       because in this case  pixcolor should not call dist.  If you do not
       have pixcolor before,  you should return 0 here. */
   if (nz < bailout) return 1; // not escaping, rare
   /*  If The Julia set is disconnected and the orbit of z comes close to
    z = 0 before escaping,  nzp will be small  */
   if (nzp < nz) return 10; // includes escaping through 0
 
   // compute estimated distance  = x 
   x = log(nz); 
   x = x*x*nz / (temp[1]*nzp); //4*square of dist/pixelwidth
   if (x < 0.04) return 1; // exterior but very close to the boundary
   if (x < 0.24) return 9; // exterior but close to the boundary
   return 10; //exterior : escaping and not close to the boundary
} //dist

Two algorithms in one loop[edit]

  • distance rendering for fractals by ińigo quilez in c++[6]


Function GiveDistance(xy2,eDx,eDy:extended):extended;
 
begin
    result:=2*log2(sqrt(xy2))*sqrt(xy2)/sqrt(sqr(eDx)+sqr(eDy));
  end;
 
//------------------------------------
 
Function PointIsInCardioid (Cx,Cy:extended):boolean;
 //Hugh Allen
 // http://homepages.borland.com/ccalvert/Contest/MarchContest/Fractal/Source/HughAllen.zip
  var DeltaX,DeltaY:extended;
      //
      PDeltyX,PDeltyY:extended;
      //
      ZFixedX,ZFixedY:extended;
 
  begin
     result:=false;
     // cardioid checkig - thx to Hugh Allen
     //sprawdzamy Czy punkt  C jest w głównej kardioidzie
     //Cardioid in squere :[-0.75,0.4] x [ -0.65,0.65]
     if InRange(Cx,-0.75,0.4)and InRange(Cy,-0.65,065) then
      begin
        // M1= all C for which Fc(z) has  attractive( = stable) fixed point
        // znajdyjemy punkt staly z: z=z*z+c
        // czyli rozwiazujemy rownanie kwadratowe
        // zmiennej zespolonej o wspolczynnikach zespolonych
        // Z*Z - Z + C = 0
        //Delta:=1-4*a*c; Delta i C sa liczbami zespolonymi
        DeltaX:=1-4*Cx;
        DeltaY:=-4*Cy;
        // Pierwiastek zespolony z delty
        CmplxSqrt(DeltaX,DeltaY,PDeltyX,PDeltyY);
        // obliczmy punkt staly jeden z dwóch, ten jest prawdopodobnie przycišgajšcy
        ZFixedX:=1.0-PDeltyX; //0.5-0.5*PDeltyX;
        ZFixedY:=PDeltyY; //-0.5*PDeltyY;
        // jesli punkt stały jest przycišgajšcy
        // to należy do M1
        If (ZfixedX*ZFixedX + ZFixedY*ZFixedY)<1.0
          then result:=true;
 
 
          // ominięcie iteracji M1 przyspiesza 3500 do 1500 msek
        end; // if InRange(Cx ...
 
 
 
 
  end;
//------------------------------------
Function PointIsInComponent (Cx,Cy:extended):boolean;
//Hugh Allen
// http://homepages.borland.com/ccalvert/Contest/MarchContest/Fractal/Source/HughAllen.zip
 
 
  var Dx:extended;
  begin
    result:=false;
    // czy punkt C nalezy do koła na lewo od kardioidy
    // circle: center = -1.0  and radius 1/4
    dx:=Cx+1.0;
    if (Dx*Dx+Cy*Cy) < 0.0625
      then result:=true;
 
 
  end;
 
//------------------------------
Procedure DrawDEM_DazibaoTrueColor; 
// draws Mandelbrot set in black and its complement in true color
// see   http://ibiblio.org/e-notes/MSet/DEstim.htm
// by  Evgeny Demidov
//
// see also
//http://www.mandelbrot-dazibao.com/Mset/Mset.htm
// translation ( with modification) of Q-Basic program:
//  http://www.mandelbrot-dazibao.com/Mset/Mdb3.bas
//
// see also my page http://republika.pl/fraktal/mset_dem.html
 
 
var iter:integer;
     iY,iX:integer;
     eCy ,eCx:extended; // C:=eCx + eCy*i
     eX,eY:extended;    // Zn:=eX+eY*i
     eTempX,eTempY:extended;
     eX2,eY2:extended;  //x2:=eX*eX;  y2:=eY*eY;
     eXY2:extended;  //  xy2:=x2+y2;
     eXY4:extended;
     eTempDx,eTempDy:extended;
     eDx,eDy:extended; // derivative
     distance:extended;
     color:TColor;
 
begin
  //compute bitmap
  for iY:= iYmin to iYMax do
    begin
      eCy:=Convert_iY_to_eY(iY);
      for iX:= iXmin to iXmax do
        begin
          eCx:=Convert_iX_to_eX(iX);
          If not PointIsInCardioid (eCx,eCy) and not PointIsInComponent(eCx,eCy)
            then
              begin
                //  Z0:=0+0*i
                eX:=0;
                eY:=0;
                eTempX:=0;
                eTempY:=0;
                //
                eX2:=0;
                eY2:=0;
                eXY2:=0;
                //
                eDx:=0;
                eDy:=0;
                eTempDx:=0;
                eTempDy:=0;
                //
                iter:=0;
                // iteration of Z ; Z= Z*z +c
                while ((iter<IterationMax) and (eXY2<=BailOut2)) do
                  begin
                    inc(iter);
                    //
                    eTempY:=2*eX*eY + eCy;
                    eTempX:=eX2-eY2 + eCx;
                    //
                    eX2:=eTempX*eTempX;
                    eY2:=eTempY*eTempY;
                    //
                    eTempDx:=1+2*(eX*eDx-eY*eDy);
                    eTempDy:=2*(eX*eDY+eY*eDx);
                    //
                    eXY2:=eX2+eY2;
                    //
                    eX:=ETempX;
                    eY:=eTempY;
                    //
                    eDx:=eTempDx;
                    eDy:=eTempDy;
                  end; // while
 
                // drawing procedure
                if (iter<IterationMax)
                then
                  begin
                    distance:= GiveDistance(eXY2,eDx,eDy);
                    color:=Rainbow(1,500,Abs(Round(100*Log10(distance)) mod 500));
                    with Bitmap1.FirstLine[iY*Bitmap1.LineLength+iX] do
                        begin
                          B := GetBValue(color);
                          G := GetGValue(color);
                          R := GetRValue(color);
                          //A := 0;
                        end; // with FirstLine[Y*LineLength+X]
 
                end // if (iter<IterationMax) then
              else  with Bitmap1.FirstLine[iY*Bitmap1.LineLength+iX] do
                        begin
                          B := 0;
                          G := 0;
                          R := 0;
                          //A := 0;
                        end;
             //--- end of drawing procedure ---
            end //  If not PointIsInCardioid ... then
          else with Bitmap1.FirstLine[iY*Bitmap1.LineLength+iX] do
                        begin
                          B := 0;
                          G := 0;
                          R := 0;
                          //A := 0;
                        end;
         //--- If not PointIsInCardioid ...
        end; // for iX
      end; // for iY
 
 
end;
//------------------------------------------------------

Optimistaion[edit]

Creating DEM images can be improved by use of :

Interior distance estimation[edit]

Estimates distance from limitly periodic point c\, ( in the interior of Mandelbrot set ) to the boundary of Mandelbrot set.


Descriptions :[10]

Pixels colored according to the estimated interior distance

Math description[edit]

The estimate is given by :


distance(c,\sigma M) = d(c, z_0, p) = \frac{1-\mid{\frac{\partial}{\partial{z}}f_c^p(z_0)}\mid^2}
  {\mid{\frac{\partial}{\partial{c}}\frac{\partial}{\partial{z}}f_c^p(z_0) +
        \frac{\partial}{\partial{z}}\frac{\partial}{\partial{z}}f_c^p(z_0)
        \frac{\frac{\partial}{\partial{c}}f_c^p(z_0)}
             {1-\frac{\partial}{\partial{z}}f_c^p(z_0)}}\mid}

where


p is the period = length of the periodic orbit

c is the point from parameter plane to be estimated

f_c(z) is the complex quadratic polynomial f_c(z)=z^2 + c

f_c^p(z_0) is the p-fold iteration of f_c(z)

z_0 is any of the p points that make the periodic orbit ( limit cycle ) \{z_0, \dots , z_{p-1} \}

4 derivatives of f_c^p evaluated at z_0, c, p :

\frac{\partial}{\partial{c}}\frac{\partial}{\partial{z}}f_c^p(z_0)

\frac{\partial}{\partial{z}}\frac{\partial}{\partial{z}}f_c^p(z_0)

\frac{\partial}{\partial{c}}f_c^p(z_0)

\frac{\partial}{\partial{z}}f_c^p(z_0)

First partial derivative with respect to z[edit]

It must be computed recursively by applying the chain rule :

Starting points : 
\begin{cases}
\frac{\partial}{\partial{z}}f_c^0(z_0) = 1 \\
f_c^0(z_0) = z_0
\end{cases}

First iteration : 
\begin{cases}
\frac{\partial}{\partial{z}}f_c^1(z_0) = 2*f_c^0(z_0)*\frac{\partial}{\partial{z}}f_c^0(z_0) =  2z_0 \\
f_c^1(z_0) = (f_c^0(z_0))^2 + c = z_0^2 + c
\end{cases}


Second iteration : 
\begin{cases}
\frac{\partial}{\partial{z}}f_c^2(z_0) =2*f_c^1(z_0)*\frac{\partial}{\partial{z}}f_c^1(z_0)= 4z_0^3+4cz_0  \\
f_c^2(z_0) = (f_c^1(z_0))^2 + c = (z_0^2 + c)^2 + c = z_0^4 + 2*c*z_0^2 + c^2 + c
\end{cases}


General rule for p iteration : 
\begin{cases}
\frac{\partial}{\partial{z}}f_c^p(z_0) = 2*f_c^{p-1}(z_0)*\frac{\partial}{\partial{z}}f_c^{p-1}(z_0)  \\
f_c^p(z_0) = f_c^{p-1}(z_0)^2 + c
\end{cases}

First partial derivative with respect to c[edit]

It must be computed recursively by applying the chain rule :

Starting points : 
\begin{cases}
\frac{\partial}{\partial{c}}f_c^p(z_0) = 0 \\
f_c^0(z_0) = z_0
\end{cases}

First iteration : 
\begin{cases}
\frac{\partial}{\partial{c}}f_c^1(z_0) = 2*f_c^0(z_0)*\frac{\partial}{\partial{c}}f_c^0(z_0) + 1 =  1 \\
f_c^1(z_0) = (f_c^0(z_0))^2 + c = z_0^2 + c
\end{cases}


Second iteration : 
\begin{cases}
\frac{\partial}{\partial{c}}f_c^2(z_0) =2*f_c^1(z_0)*\frac{\partial}{\partial{c}}f_c^1(z_0) + 1 = 2z_0^2 +2c + 1 \\
f_c^2(z_0) = (f_c^1(z_0))^2 + c = (z_0^2 + c)^2 + c = z_0^4 + 2*c*z_0^2 + c^2 + c
\end{cases}


General rule for p iteration : 
\begin{cases}
\frac{\partial}{\partial{c}}f_c^p(z_0) = 2*f_c^{p-1}(z_0)*\frac{\partial}{\partial{c}}f_c^p(z_0) + 1  \\
f_c^p(z_0) = f_c^{p-1}(z_0)^2 + c
\end{cases}

Second order partial derivative with respect to z[edit]

It must be computed :

  • together with :f_c^n and \frac{\partial}{\partial{z}}f_c^p(z_0)
  • recursively by applying the chain rule

Starting points : 
\begin{cases}
\frac{\partial}{\partial{z}}f_c^0(z_0) = 1 \\
\frac{\partial}{\partial{z}}\frac{\partial}{\partial{z}}f_c^p(z_0) = 0 \\
f_c^0(z_0) = z_0
\end{cases}

First iteration : 
\begin{cases}
\frac{\partial}{\partial{z}}f_c^1(z_0) = 2*f_c^0(z_0)*\frac{\partial}{\partial{z}}f_c^0(z_0) =  2z_0 \\
\frac{\partial}{\partial{z}}\frac{\partial}{\partial{z}}f_c^p(z_0) = 2\{ (\frac{\partial}{\partial{z}}f_c^0(z_0))^2 + f_c^0(z_0)*\frac{\partial}{\partial{z}}\frac{\partial}{\partial{z}}f_c^0(z_0)\} = 2\\
f_c^1(z_0) = (f_c^0(z_0))^2 + c = z_0^2 + c
\end{cases}


Second iteration : 
\begin{cases}
\frac{\partial}{\partial{z}}f_c^2(z_0) =2*f_c^1(z_0)*\frac{\partial}{\partial{z}}f_c^1(z_0)= 4z_0^3+4cz_0  \\
\frac{\partial}{\partial{z}}\frac{\partial}{\partial{z}}f_c^p(z_0) = \\
f_c^2(z_0) = (f_c^1(z_0))^2 + c = (z_0^2 + c)^2 + c = z_0^4 + 2*c*z_0^2 + c^2 + c
\end{cases}


General rule for p iteration : 
\begin{cases}
\frac{\partial}{\partial{z}}f_c^p(z_0) = 2*f_c^{p-1}(z_0)*\frac{\partial}{\partial{z}}f_c^{p-1}(z_0)  \\
\frac{\partial}{\partial{z}}\frac{\partial}{\partial{z}}f_c^p(z_0) =  2\{ (\frac{\partial}{\partial{z}}f_c^{p-1}(z_0))^2 + f_c^{p-1}(z_0)*\frac{\partial}{\partial{z}}\frac{\partial}{\partial{z}}f_c^{p-1}(z_0)\}\\
f_c^p(z_0) = f_c^{p-1}(z_0)^2 + c
\end{cases}

Second order mixed partial derivative[edit]

\frac{\partial}{\partial{c}}\frac{\partial}{\partial{z}}f_c^p(z_0)

Algorithm[edit]

Steps[edit]

  • choose c
  • check if critical point for given c is periodic ( interior ) or not ( exterior or near boundary or to big period )
    • if not periodic do not use this algorithm ( use exterior version)
    • if periodic compute period and periodic point z_0
  • using  c, z_0, p compute distance using p iteration

Computing period and periodic point[edit]

Computing distance[edit]

Code[edit]

  • Java code by Albert Lobo Cusidó[11]
  • Haskell code by Claude[12] and image [13]

Problems[edit]

There are two practical problems with the interior distance estimate: first, we need to find z_0 precisely, and second, we need to find p precisely. The problem with z_0 is that the convergence to z_0 by iterating P_c(z) requires, theoretically, an infinite number of operations. The problem with period is that, sometimes, due to rounding errors, a period is falsely identified to be an integer multiple of the real period (e.g., a period of 86 is detected, while the real period is only 43=86/2). In such case, the distance is overestimated, i.e., the reported radius could contain points outside the Mandelbrot set.

Optimisation[edit]

Analogous to the exterior case, once b is found, we know that all points within the distance of b/4 from c are inside the Mandelbrot set.

References[edit]

  1. fractalforums discussion : How are B&W images from "Beauty of Fractals" created?
  2. Milnor (in Dynamics in One Complex Variable, appendix G)
  3. 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
  4. distance rendering for fractals by ińigo quilez
  5. Distance Estimator by Robert P. Munafo
  6. distance rendering for fractals by ińigo quilez
  7. distance estimation by Claude Heiland-Allen
  8. distance rendering for fractals by ińigo quilez
  9. fractalforums gallery by Pauldelbrot
  10. Interior and exterior distance bounds for the Mandelbrot by Albert Lobo Cusidó
  11. Interior and exterior distance bounds for the Mandelbrot by Albert Lobo Cusidó
  12. Interior coordinates in the Mandelbrot set by Claude
  13. Fractal forum : Coloring points inside the Mandelbrot set