Fractals/Iterations in the complex plane/Mandelbrot set/boundary

=Intro= Boundary of Mandelbrot set consist of :
 * boundaries of primitive and satellite hyperbolic components of Mandelbrot set including points :
 * Parabolic (including 1/4 and primitive roots which are landing points for 2 parameter rays with rational external angles = biaccesible ).
 * Siegel ( a unique parameter ray landing with irrational external angle)
 * Cremer ( a unique parameter ray landing with irrational external angle)
 * Boundary of M without boundaries of hyperbolic components with points
 * non-renormalizable (Misiurewicz with rational external angle and other).
 * finitely renormalizable (Misiurewicz and other).
 * infinitely renormalizble (Feigenbaum and other). Angle in turns of external rays landing on the Feigenbaum point are irrational numbers
 * boundaries of non hyperbolic components ( we believe they do not exist but we cannot prove it ) Boundaries of non-hyperbolic components would be infinitely renormalizable as well.



=compute =

compute t from c
Function describing relation between parameter c and internal angle t :

$$c = c(t) = \frac\lambda2\left(1-\frac\lambda2\right)= \frac{e^{\Pi*t*i}}{2} - \frac{e^{2*\Pi*t*i}}{4} $$ It is used for computing c point of boundary of main cardioid

To compute t from c one can use Maxima CAS: (%i1) eq1:c = exp(%pi*t*%i)/2 - exp(2*%pi*t*%i)/4;

%i %pi t    2 %i %pi t                             %e           %e (%o1)                   c = -- - 2            4 (%i2) solve(eq1,t); %i log(1 - sqrt(1 - 4 c))       %i log(sqrt(1 - 4 c) + 1) (%o2) [t = - -, t = - -] %pi                             %pi

So :

f1(c):=float(cabs( - %i* log(1 - sqrt(1 - 4* c))/%pi)); f2(c):=float(cabs( - %i* log(1 + sqrt(1 - 4* c))/%pi));

=Drawing boundaries=

Methods used to draw boundary of Mandelbrot set :
 * whole boundary
 * Distance Estimation Method for Mandelbrots set = DEM/M ( exterior, interior or both)
 * Boundary scanning method (BSM) : detecting the edges after Boolean escape time
 * boundary tracing method
 * contour integration by Robert Davies
 * Jungreis function $$\Psi_M\,$$
 * implicit polynomial curves that converge to the boundary of the Mandelbrot set ( lemniscates)
 * drawing boundaries of components of Mandelbrot set
 * Newton method

Jungreis function


Description :
 * Drawing Mc by Jungreis Algorithm
 * Maxima CAS src code
 * On Benford's Law and the Coefficients of the Riemann Mapping Function for the Exterior of the Mandelbrot Set 8 Jun 2022 ·  Filippo Beretta, Jesse Dimino, Weike Fang, Thomas C. Martinez, Steven J. Miller, Daniel Stoll and code

Python code 

How to draw boundaries of hyperbolic components
Methods of drawing boundaries:
 * solving boundary equations :
 * using method by Brown,Stephenson,Jung . It  works only up to period 7-8
 * using resultants
 * parametrisation of boundary with Newton method near centers of components . This methods needs centers, so it has some limitations,


 * Boundary scanning : detecting the edges after detecting period by checking every pixel. This method is slow but allows zooming. Draws "all" components
 * Boundary tracing for given c value. Draws single component.
 * Fake Mandelbrot set by Anne M. Burns : draws main cardioid and all its descendants. Do not draw mini Mandelbrot sets.

"... to draw the boundaries of hyperbolic components using Newton's method. That is, take a point in the hyperbolic component that you are interested in (where there is an attracting cycle), and then find a curve along which the modulus of the multiplier tends to one. Then you will have found an indifferent parameter. Now you can similarly change the argument of the multiplier, again using Newton's method, and trace this curve. Some care is required near "cusps". " Lasse Rempe-Gillen

solving boundary equations
/* maxima CAS*/ kill(all);

f(z):=z^2+c;

d:diff(f(z),z,1);

r:1; /* only boundary */

/* system of equations */ e2:d = r*exp(2*%pi*%i*t); e1:f(z) = z;

ss: solve([e1,e2],[z,c]); s: ss[1]; s:s[2]; s:rhs(s);

/* change form */ x:realpart(s); y:imagpart(s);

load(draw); draw2d( line_type = solid,  nticks= 500,  parametric(x,y, t,0,2*%pi) );

System of 2 equations defining boundaries of period n hyperbolic components

 * first defines periodic point,
 * second defines indifferent orbit ( multiplier of periodic point is equal to one ).

$$\begin{cases} F(n,z,c) = z \\ abs(\lambda (z) )=1 \end{cases}$$

Because stability index $$abs(\lambda)\, $$ is equal to radius of point of unit circle $$abs(w)\,$$:

$$abs(\lambda) = abs(w)\,$$

so one can change second equation to form :

$$\lambda = w \,$$

It gives system of equations :

$$\begin{cases} F(n,z,c) = z \\ \lambda = w \end{cases}$$

It can be used for :
 * drawing components ( boundaries, internal rays )
 * finding indifferent parameters ( parabolic or for Siegel discs )

Above system of 2 equations has 3 variables : $$z,c,w \,$$ ( $$n\,$$ is constant and multiplier $$\lambda\,$$ is a function of $$z,c \,$$). One have to remove 1 variable to be able to solve it.

Boundaries are closed curves : cardioids or circles. One can parametrize points of boundaries with angle $$t\,$$ ( here measured in turns from 0 to 1 ).

After evaluation of $$w = l(t) \,$$ one can put it into above system, and get a system of 2 equations with 2 variables $$z,c\,$$.

Now it can be solved

For periods:


 * 1-3 it can be solved with symbolical methods and give implicit ( boundary equation) $$b_p(w,c)=0 \,$$ and explicit function (inverse multiplier map) : $$c=\gamma_p(w)\,$$
 * 1-2 it is easy to solve
 * 3 it can be solve using "elementary algebra" ( Stephenson )
 * >3 it can't be reduced to explicitly function but
 * can be reduced to implicit solution ( boundary equation) $$b_p(w,c)=0 \,$$ and solved numerically
 * can be solved numerically using Newton method for solving system of nonlinear equations

Examples
 * . An introduction to the Mandelbrot set II May 21, 2021, by Juan Rivera-Letelier Slides, Sage notebook,

Solving system of equation for period 1
Here is Maxima code :

(%i4) p:1; (%o4) 1 (%i5) e1:F(p,z,c)=z; (%o5) z^2+c=z (%i6) e2:m(p)=w; (%o6) 2*z=w (%i8) s:eliminate ([e1,e2], [z]); (%o8) [w^2-2*w+4*c] (%i12) s:solve([s[1]], [c]); (%o12) [c=-(w^2-2*w)/4] (%i13) define (m1(w),rhs(s[1])); (%o13) m1(w):=-(w^2-2*w)/4

where w:exp(2*%pi*%i*t); and (%i15) display2d:false; (%o15) false (%i16) 2*w; (%o16) 2*%e^(2*%i*%pi*t) (%i17) w*w; (%o17) %e^(4*%i*%pi*t)

Second equation contains only one variable, one can eliminate this variable. Because boundary equation is simple so it is easy to get explicit solution

m1(w):=-(w^2-2*w)/4

So

m1(t) := block([w], w:exp(2*%pi*%i*t), return ((2*w-w^2)/4)); g(t) := float(rectform(t));

Math equation:

$$c_t = \frac{2 e^{2\pi it} - e^{4\pi it }}{4}$$

Solving system of equation for period 2
Here is Maxima code using to_poly_solve package by Barton Willis:

(%i4) p:2; (%o4) 2 (%i5) e1:F(p,z,c)=z; (%o5) (z^2+c)^2+c=z (%i6) e2:m(p)=w; (%o6) 4*z*(z^2+c)=w (%i7) e1:F(p,z,c)=z; (%o7) (z^2+c)^2+c=z (%i10) load(topoly_solver); to_poly_solve([e1, e2], [z, c]); (%o10) C:/PROGRA~1/MAXIMA~1.1/share/maxima/5.16.1/share/contrib/topoly_solver.mac (%o11) z=sqrt(w)/2,c=-(w-2*sqrt(w))/4],[z=-sqrt(w)/2,c=-(w+2*sqrt(w))/4],[z=(sqrt(1-w)-1)/2,c=(w-4)/4],[z=-(sqrt(1-w)+1)/2,c=(w-4)/4 (%i12) s:to_poly_solve([e1, e2], [z, c]); (%o12) z=sqrt(w)/2,c=-(w-2*sqrt(w))/4],[z=-sqrt(w)/2,c=-(w+2*sqrt(w))/4],[z=(sqrt(1-w)-1)/2,c=(w-4)/4],[z=-(sqrt(1-w)+1)/2,c=(w-4)/4 (%i14) rhs(s[4][2]); (%o14) (w-4)/4 (%i16) define (m2 (w),rhs(s[4][2])); (%o16) m2(w):=(w-4)/4

explicit solution :

m2(w):=(w-4)/4

Solving system of equation for period 3
For period 3 ( and higher) previous method give no results (Maxima code) :

(%i14) p:3; e1:z=F(p,z,c); e2:m(p)=w; load(topoly_solver); to_poly_solve([e1, e2], [z, c]); (%o14) 3 (%o15) z=((z^2+c)^2+c)^2+c (%o16) 8*z*(z^2+c)*((z^2+c)^2+c)=w (%i17) (%o17) C:/PROGRA~1/MAXIMA~1.1/share/maxima/5.16.1/share/contrib/topoly_solver.mac `algsys' cannot solve - system too complicated. #0: to_poly_solve(e=[z = ((z^2+c)^2+c)^2+c,8*z*(z^2+c)*((z^2+c)^2+c) = w],vars=[z,c]) -- an error. To debug this try debugmode(true);

I use code by Robert P. Munafo which is based on paper of Wolf Jung.

One can approximate period 3 components with equations :

(%i1) z:x+y*%i; (%o1)                             %i y + x (%i2) w:asinh(z); (%o2)                          asinh(%i y + x) (%i3) realpart(w); (%o3) 2   2          2    2     2      2  2 1/4     atan2(2 x y, - y  + x  + 1)      2 log((((- y + x  + 1)  + 4 x  y )    sin(---) + y)                                                      2                                                        2    2         2    2     2      2  2 1/4     atan2(2 x y, - y  + x  + 1)      2 + (((- y  + x  + 1)  + 4 x  y )    cos(---) + x) )/2 2 (%i4) imagpart(w); 2   2                 2    2     2      2  2 1/4     atan2(2 x y, - y  + x  + 1) (%o4) atan2(((- y + x  + 1)  + 4 x  y )    sin(---)                                                             2                                                              2    2               2    2     2      2  2 1/4     atan2(2 x y, - y  + x  + 1)     + y, ((- y  + x  + 1)  + 4 x  y )    cos(---) + x)

Boundary equation
The result of solving above system with respect to $$c\,$$ is boundary equation,


 * $$b_p(w,c)=0 \,$$

where $$b_p(w,c)\,$$ is boundary polynomial.

It defines exact coordinates of hyperbolic componenets for given period $$p \,$$.

It is implicit equation.

Period 1

One can easly compute boundary point c

$$ c = c_x + c_y*i $$

of period 1 hyperbolic component ( main cardioid) for given internal angle ( rotation number) t using this code by Wolf Jung

Period 2

Solving boundary equations
Solving boundary equations for various angles gives list of boundary points.

You can use Newton's method in 2 complex variables to find points in a particular component with a given multiplier value:
 * m_d_interior from mandelbrot-numerics library

Use it something like:

double _Complex z = 0; double _Complex c = /* nucleus of component */; int period = /* period of component */; for (int t = 0; t < 360; ++t) { m_d_interior(&z, &c, z, c, cexp(2 * M_PI * I * (t + 0.5) / 360), period, 16); /* plot c */ }

= size =

size of the component

 * estimated size of componnet
 * on_the_precision_required_for_size_estimates by Claude Heiland-Allen
 * deriving_the_size_estimate by Claude Heiland-Allen
 * the N-2 rule
 * The approximation ( above ) was proved by Guckenheimer and McGehee : Guckenheimer, J., and McGehee, R., "a proof of the mandelbrot n2 conjecture," institut mittag-leffler, report 15, 1984.
 * Bifurcation points and the "n-squared" approximation
 * atom domain size

=distinguish cardioids from pseudocircles= method of distinguish cardioids from pseudocircles is described in : Universal Mandelbrot Set by A. Dolotin. The relevant part is section 5, in particular equation 5.8. In the paper $$\gamma$$ is defined implicitly in equation 5.1,

$$\gamma = \frac{\partial}{\partial z} f(z, c)$$

Equation 5.8 then becomes: $$ \epsilon = - \frac{1}{ \frac{\partial}{\partial c}F \frac{\partial}{\partial z}F} \left(\frac{\frac{\partial}{\partial c}\frac{\partial}{\partial c}F}{2 \frac{\partial}{\partial c}F} + \frac{\frac{\partial}{\partial c}\frac{\partial}{\partial z}F}{\frac{\partial}{\partial z}F}\right) $$

When epsilon is:
 * near 0 then you have a cardioid
 * near 1 then you have a circle

Code: =distorsion and eccentricity=
 * c by Claude
 * In general, the eccentricity grows with increasing period, but not linearly, not even monotone. It looks from the first 40 periods as if they approach perhaps already a limit value (image 2).
 * fractalforums.com : how-distorted-can-a-minibrot-be ?

=degree n =
 * The boundary of the central hyperbolic component of Md is an epicycloid with d − 1 cusps
 * boundary of period 2 component
 * Locating and counting bifurcation points of satellite components from the main component in the degree- n bifurcation set
 * Accurate Computation of Periodic Regions’ Centers in the General M-Set

=References=
 * matheplanet article
 * fractalforums : how-would-you-define-a-path-along-mandelbrot-border