Fractal Software > Programming

exterior distance estimation for logistic map

(1/7) > >>

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

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$

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.

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

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