Fractal Software > Programming

 exterior distance estimation for logistic map

(1/7) > >>

Adam Majewski:
Hi

I try to use distance estimation for Mandelbrot set of the logistic map (DEM/M )

\[ f(z,m) = m*z*(1-z) \]

Iterated function :

\[ f^0(z_0,m) = z_0 \]
\[ f^{n+1}(z_0,m) = f(f^n(z_0, m), m) \]
\[ f(z_0,m) = f^1(z_0, m) \]



For DEM/M I need to compute first derivative of f wrt to m.

For simple function \( f \) it seems easy :

\[ \dfrac{\partial}{\partial m} f(z, m) = z - z^2 \]


For iterated function there are rules:

* Iterated function \( f^{n+1} (z, m) = f (f^n (z, m), m) \) is a function composition
* Applay chain rule for composite functions, so if  \( f(x) = h(g(x)) \) then \( f'(x) = h'(g(x)) \cdot g'(x) \)
* \( z_0 = 0.5 \) is a constant value
* m is a variable here

I use arbitrary names

* \( A_n = f^n(z_0, m) \)
* \( D_n = \dfrac{\partial}{\partial m} f^n (z_0, m) \)
I compute initial values first

\[ A_0 \quad \xrightarrow{definition\ of\ A}\  f^0 (z_0,m) \quad \xrightarrow{definition\ of\ f}\ \quad z_0  \]

\[ D_0 \quad \xrightarrow{definition\ of\ D } \dfrac{\partial}{\partial m} f^0 (z_0, m)  \xrightarrow{ definition\ of\  f } \dfrac{\partial}{\partial m} z_0 \xrightarrow{ derivative\ of\ constant\ value} 0 \]

next values :

\[ A_{n+1}  \quad \xrightarrow{definition\ of\ A} f^{n+1} (z_0, m)  \quad \xrightarrow{definition\ of\ f} f (f^n (z_0, m), m)  \quad \xrightarrow{definition\ of\ A} f (A_n , m)  \quad \xrightarrow{definition\ of\ f} m*A_n  - m*A_n^2  \]

\[ D_{n+1}  \quad \xrightarrow{definition\ of\ D} \dfrac{\partial}{\partial m} A_{n+1}  \quad \xrightarrow{definition\ of\ A} \dfrac{\partial}{\partial m}  f(f^n (z_0, m), m) \quad \xrightarrow{chain\ rule} (m*A_n - m*A_n^2)* D_n   \]

Am I right ?
TIA


Linkback: https://fractalforums.org/index.php?topic=4458.0

marcm200:
I think the last step for D_(n+1) should be
\[
D_{n+1}=\frac{d}{dm} A_{n+1} = \frac{d}{dm} f(A_n)=f'(A_n)\cdot \frac{d}{dm} A_n = (A_n-A_n^2)\cdot D_n
 \]



Adam Majewski:
but if \( D_0= 0 \) then \( D_n  \) will allways be zero.

???

marcm200:
Why go down to D0? In iterations you never do the 0-step iteration but always at least 1 to get fix-points. So the induction start for me would be f1, D1 and A1.

Secondly, f[0-fold iterated] is an ambiguous definition. It is the identity function. But what is the identity function of a two-variate function f(z,m)? If you're working with d/dm it is imo f0(z,m)=m , and if you're working with z it is f0(z,m)=z, so the derivative is 1.

Adam Majewski:
Ok.

I have used dem for z^2+c and m*z(1-z) . Here are 2 images . It seems that it works for c but not for m = does not find boundary

Here is the diff:

--- Code: ---73,74c73,74
< const int iXmax = 800;
< const int iYmax = 801;
---
> const int iXmax = 1600;
> const int iYmax = 800;
81,84c81,84
< const double CxMin=-2.25;
< const double CxMax=0.75;
< const double CyMin=1.5;
< const double CyMax= -1.5;
---
> const double CxMin= -1.6;
> const double CxMax=  3.6;
> const double CyMin= -1.35;
> const double CyMax=  1.35;
93c93
< char *filename="demmm_i090.ppm"; //
---
> char *filename="demm_log.ppm"; //
139c139
<   double complex Z= 0.0; // initial value for iteration Z0
---
>   double complex Z= 0.5; // initial value for iteration Z0
142c142
<   double complex dC = 0; // derivative
---
>   double complex dC = 1.0; // derivative
150,151c150,151
<       dC = 2 * Z * dC + 1;
<       Z=Z*Z+C; // https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/qpolynomials
---
>       dC =  Z *(1-Z)* dC ;
>       Z=Z*(1 -Z)*C; // https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/qpolynomials
243c243,244
<   
---
>   printf("PixelWidth = %.16f \n", PixelWidth );
>    printf("PixelHeight = %.16f \n", PixelHeight );

--- End code ---


and the code for m

--- Code: ---
/*
   c program: console, 1-file
   d.c
   
   algorithms:
   - escape time
   - DEM/M = distance estimation, Milnor version
   
   
   "we have an estimate dn for the distance from z0 to M,
   where n is the first index for which |zn|>ER.
   ...
   we compared this estimate to a fixed size s and set the pixel
   * black if dn<s
   * white otherwise" A Cheritat
   
   
 
   
   --------------------------------
   1. draws Mandelbrot set for complex quadratic polynomial
   Fc(z)=z*z +c
   using samm = Stripe Average Method/Coloring 
   -------------------------------         
   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)
   -----
   it is example  for :
   https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set
 
   -------------
   compile :

 
 
   gcc d.c -lm -Wall
 
 
   ./a.out
   
   convert samm.ppm -resize 800x800 samm.png

   -------- git --------
   
   
   cd existing_folder
   git init
   git remote add origin git@gitlab.com:adammajewski/mandelbrot_wiki_ACh.git
   git add samm.c
   git commit -m "samm + linear interpolation "
   git push -u origin master



 
 
 
*/
#include <stdio.h>
#include <math.h>
#include <complex.h> // https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c

 


 
/* screen ( integer) coordinate */
int iX,iY;
const int iXmax = 1600;
const int iYmax = 800;

/* world ( double) coordinate = parameter plane*/
double Cx,Cy;

// c = -0.355298979690605  +1.178016112830515 i    period = 0
// c = 0.213487557597331  +0.604701020309395 i    period = 0
const double CxMin= -1.6;
const double CxMax=  3.6;
const double CyMin= -1.35;
const double CyMax=  1.35;
/* */
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="demm_log.ppm"; //
char *comment="# ";/* comment should start with # */
       
static unsigned char color[3]; // 24-bit rgb color


const int IterationMax=2500; //  N in wiki


/* bail-out value for the bailout test for exaping points
   radius of circle centered ad the origin, exterior of such circle is a target set  */
const double EscapeRadius=1000; // = ER big !!!!

double thickness = 0.90;
double MinBoundaryWidth;       
       
       
       
double complex give_c(int iX, int iY){
  double Cx,Cy;
 
  Cy=CyMax - iY*PixelHeight; // reverse Y axis
  //if (fabs(Cy)<PixelHeight) Cy = 0.0;
  Cx=CxMin + iX*PixelWidth;
   
  return Cx+Cy*I;
 
 
}
 




/*
  input :
  - complex number
  - intege
  output =  estimator
   
*/
double Give_D(double complex C , int iMax)
{
  int i=0; // iteration
   
   
  double complex Z= 0.5; // initial value for iteration Z0
  double R; // =radius = cabs(Z)
  double D;
  double complex dC = 1.0; // derivative
  double de; // = 2 * z * log(cabs(z)) / dc;
   
   
  // iteration = computing the orbit
  for(i=0;i<iMax;i++)
    {
   
      dC =  Z *(1-Z)* dC ;
      Z=Z*(1 -Z)*C; // https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/qpolynomials
     
           
      R = cabs(Z);
      if(R > EscapeRadius) break; // exterior of M set
   
     
    } // for(i=0
   
   
  if (i == iMax) D = -1.0; // interior
    else { // exterior
      de = 2.0 * R * log(R) / cabs(dC) ; //
   
      if (de < MinBoundaryWidth) D = FP_ZERO; //  boundary
         else  D = 1.0; // exterior
    }
   
  return D; 
}
 

/*
 
 input = complex number
 output
  - color array as Formal parameters
  - int = return code
*/
int compute_color(complex double c, unsigned char color[3]){
 
  double d ;
  unsigned char b;
   
  // compute
  d = Give_D( c, IterationMax);
     
  //
  if (d < 0.0)
     /*  interior of Mandelbrot set = inside_color =  */
      b = 191;                         
   
  else // exterior and boundary
     {     
       if (d == FP_ZERO) b = 255; // boundary     
          else b = 0; // exterior
     };
   
 // gray gradient
 color[0]= b;  /* Red*/
 color[1]= b;  /* Green */
 color[2]= b;  /* Blue */
 
 
  return 0;
}
 
 
 
void setup(){
 
  //
  PixelWidth=(CxMax-CxMin)/iXmax;
  PixelHeight=(CyMax-CyMin)/iYmax;
  //
  MinBoundaryWidth = thickness*PixelWidth;
 
       
         
  /*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);
 
}
 



void info(){

  double distortion;
  // widt/height
  double PixelsAspectRatio = (double)iXmax/iYmax;  // https://en.wikipedia.org/wiki/Aspect_ratio_(image)
  double WorldAspectRatio = (CxMax-CxMin)/(CyMax-CyMin);
  printf("PixelsAspectRatio = %.16f \n", PixelsAspectRatio );
  printf("WorldAspectRatio = %.16f \n", WorldAspectRatio );
  distortion = PixelsAspectRatio - WorldAspectRatio;
  printf("distortion = %.16f ( it should be zero !)\n", distortion );
  printf("\n");
  printf("bailut value = Escape Radius = %.0f \n", EscapeRadius);
  printf("IterationMax = %d \n", IterationMax);
  printf("PixelWidth = %.16f \n", PixelWidth );
   printf("PixelHeight = %.16f \n", PixelHeight );

  // file 
  printf("file %s saved. It is called MilMand_1000iter in  A Cheritat wiki\n", filename);

}
 



 
void close(){
 
  fclose(fp);
  info();

}
 
 
 
 
 
// ************************************* main *************************
int main()
{
       
  complex double c;
       
       
 
  setup();     
       
       
  printf(" render = compute and write image data bytes to the file \n");
 //#pragma omp parallel for schedule(static) private(iY, iX, c, color)
  for(iY=0;iY<iYmax;iY++)
    for(iX=0;iX<iXmax;iX++)
      { // compute pixel coordinate       
c = give_c(iX, iY); 
/* compute  pixel color (24 bit = 3 bytes) */
compute_color(c,color);         
/*write color to the file*/
fwrite(color,1,3,fp);
      }
       
 
 
  close();
 
         
  return 0;
}


--- End code ---

Navigation

[0] Message Index

[#] Next page

Go to full version
Powered by SMFPacks Rates (Facepunch) Mod
Powered by SMFPacks SEO Pro Mod | Sitemap
Powered by Advanced Topic Prefix Pro
Powered by SMFPacks Advanced Attachments Uploader Mod