# Fractals/Iterations in the complex plane/periodic points

## period

Here period means period of iterated function (map).

Point z has period p under f if :

$z:\ f^{p}(z)=z$ or in other notation:

$z_{p}=z_{0}$ The smallest positive integer p satisfying the above is called the prime period or least period of the point z.

## equation and polynomial F

So equation for periodic points:

$f^{p}(z)=z$ or in other notation

$F^{p}(z)=f^{p}(z)-z=0$ ## degree

Degree d of polynomial F :

$d=2^{p-1}$ # How to find periodic points for given period?

## Visual check

Visual check of critical orbit helps to:

• understand dynamics
• find periodic points

## Symbolic methods

 f(z,c):=z*z+c;

 fn(p, z, c) :=
if p=0 then z
elseif p=1 then f(z,c)
else f(fn(p-1, z, c),c);


### Standard polynomial $F_{p}\,$ which roots are periodic z-points of period p and its divisors

• math notation : $\ F_{p}(z,c)=f_{c}^{(p)}(z)-z$ • Maxima CAS function :
F(p, z, c) := fn(p, z, c) - z ;


### Function for computing reduced polynomial $G_{p}\,$ which roots are periodic z-points of period p without its divisors

• math definition : $G_{p}(z,c)={\frac {F_{p}(z,c)}{\prod _{m|p,m • Maxima function:
G[p,z,c]:=
block(
[f:divisors(p),
t:1], /* t is temporary variable = product of Gn for (divisors of p) other than p */
f:delete(p,f), /* delete p from list of divisors */
if p=1
then return(F(p,z,c)),
for i in f do
t:t*G[i,z,c],
g: F(p,z,c)/t,
return(ratsimp(g))
)$ ##### Finding equations for periodic points ###### Standard equations Equation for periodic point for period $p\,$ is : $\ F_{p}(z,c)=0\,$ so in Maxima : for i:1 thru 4 step 1 do disp(i,expand(Fn(i,z,c)=0)); 1 z^2-z+c=0 2 z^4+2*c*z^2-z+c^2+c=0 3 z^8+4*c*z^6+6*c^2*z^4+2*c*z^4+4*c^3*z^2+4*c^2*z^2-z+c^4+2*c^3+c^2+c=0 4 z^16+8*c*z^14+28*c^2*z^12+4*c*z^12+56*c^3*z^10+24*c^2*z^10+70*c^4*z^8+60*c^3*z^8+6*c^2*z^8+2*c*z^8+56*c^5*z^6+80*c^4*z^6+ 24*c^3*z^6+8*c^2*z^6+28*c^6*z^4+60*c^5*z^4+36*c^4*z^4+16*c^3*z^4+4*c^2*z^4+8*c^7*z^2+24*c^6*z^2+24*c^5*z^2+16*c^4*z^2+ 8*c^3*z^2- z+c^8+4*c^7+6*c^6+6*c^5+5*c^4+2*c^3+c^2+c=0  Degree of standard equations is $2^{p}\,$ These equations give periodic points for period p and its divisors. It is because these polynomials are product of lower period polynomials $G_{p}\,$ : for i:1 thru 4 step 1 do disp(i,factor(Fz(i,z,c))); 1 z^2-z+c 2 (z^2-z+c)*(z^2+z+c+1) 3 (z^2-z+c)*(z^6+z^5+3*c*z^4+z^4+2*c*z^3+z^3+3*c^2*z^2+3*c*z^2+z^2+c^2*z+2*c*z+z+c^3+2*c^2+c+1) 4 (z^2-z+c)*(z^2+z+c+1)(z^12+6*c*z^10+z^9+15*c^2*z^8+3*c*z^8+4*c*z^7+20*c^3*z^6+12*c^2*z^6+z^6+6*c^2*z^5+2*c*z^5+15*c^4*z^4+18*c^3*z^4+3* c^2*z^4+4*c*z^4+4*c^3*z^3+4*c^2*z^3+z^3+6*c^5*z^2+12*c^4*z^2+6*c^3*z^2+5*c^2*z^2+c*z^2+c^4*z+2*c^3*z+c^2*z+2*c*z+c^6+ 3*c^5+3*c^4+3*c^3+2*c^2+1)  so $\ F_{1}(z,c)=z^{2}-z+c=G_{1}\,$ $\ F_{2}(z,c)=(z^{2}-z+c)*(z^{2}+z+c+1)=G_{1}*G_{2}\,$ $\ F_{3}(z,c)=(z^{2}-z+c)*(z^{6}+z^{5}+3*c*z^{4}+z^{4}+2*c*z^{3}+z^{3}+3*c^{2}*z^{2}+3*c*z^{2}+z^{2}+c^{2}*z+2*c*z+z+c^{3}+2*c^{2}+c+1)=G_{1}*G_{3}\,$ $\ F_{4}(z,c)=G_{1}*G_{2}*G_{4}\,$ $\ F_{5}(z,c)=G_{1}*G_{5}\,$ Standard equations can be reduced to equations $G_{p}\,$ giving periodic points for period p without its divisors. See also : prime factors in wikipedia ###### Reduced equations $G_{p}(z,c)=0\,$ so in Maxima gives : for i:1 thru 5 step 1 do disp(i,expand(G[i,z,c]=0)); 1 z^2-z+c=0 2 z^2+z+c+1=0 3 z^6+z^5+3*c*z^4+z^4+2*c*z^3+z^3+3*c^2*z^2+3*c*z^2+z^2+c^2*z+2*c*z+z+c^3+2*c^2+c+1=0 4 z^12+6*c*z^10+z^9+15*c^2*z^8+3*c*z^8+4*c*z^7+20*c^3*z^6+12*c^2*z^6+z^6+6*c^2*z^5+2*c*z^5+15*c^4*z^4+18*c^3*z^4+ 3*c^2*z^4+4*c*z^4+4*c^3*z^3+4*c^2*z^3+z^3+6*c^5*z^2+12*c^4*z^2+ 6*c^3*z^2+5*c^2*z^2+c*z^2+c^4*z+2*c^3*z+c^2*z+2*c*z+c^6+3*c^5+3*c^4+3*c^3+2*c^2+1=0  ##### Computing periodic points Definition in wikipedia  ###### Fixed points ( period = 1 ) Equation defining fixed points :  z^2-z+c = 0  in Maxima : Define c value (%i1) c:%i; (%o1) %i  Compute fixed points (%i2) p:float(rectform(solve([z*z+c=z],[z]))); (%o2) [z=0.62481053384383*%i-0.30024259022012,z=1.30024259022012-0.62481053384383*%i]  Find which is repelling : abs(2*rhs(p));  Show that sum of fixed points is 1 (%i10) p:solve([z*z+c=z], [z]); (%o10) [z=-(sqrt(1-4*%i)-1)/2,z=(sqrt(1-4*%i)+1)/2] (%i14) s:radcan(rhs(p+p)); (%o14) 1  Draw points : (%i15) xx:makelist (realpart(rhs(p[i])), i, 1, length(p)); (%o15) [-0.30024259022012,1.30024259022012] (%i16) yy:makelist (imagpart(rhs(p[i])), i, 1, length(p)); (%o16) [0.62481053384383,-0.62481053384383] (%i18) plot2d([discrete,xx,yy], [style, points]);  One can add point z=1/2 to the lists: (%i31) xx:cons(1/2,xx); (%o31) [1/2,-0.30024259022012,1.30024259022012] (%i34) yy:cons(0,yy); (%o34) [0,0.62481053384383,-0.62481053384383] (%i35) plot2d([discrete,xx,yy], [style, points]);  In C : double complex GiveAlfaFixedPoint(double complex c) { // Equation defining fixed points : z^2-z+c = 0 // a = 1, b = -1 c = c double dx, dy; //The discriminant is the d=b^2- 4ac = dx+dy*i double r; // r(d)=sqrt(dx^2 + dy^2) double sx, sy; // s = sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx + sy*i double ax, ay; // alfa = ax+ay*I // d=1-4c = dx+dy*i dx = 1 - 4*creal(c); dy = -4 * cimag(c); // r(d)=sqrt(dx^2 + dy^2) r = sqrt(dx*dx + dy*dy); //sqrt(d) = s =sx +sy*i sx = sqrt((r+dx)/2); sy = sqrt((r-dx)/2); // alfa = ax +ay*i = (1-sqrt(d))/2 = (1-sx + sy*i)/2 ax = 0.5 - sx/2.0; ay = sy/2.0; return ax+ay*I; }  ###### period 2 points Using non-reduced equation in Maxima CAS : (%i2) solve([(z*z+c)^2+c=z], [z]); (%o2) [z=-(sqrt(-4*c-3)+1)/2,z=(sqrt(-4*c-3)-1)/2,z=-(sqrt(1-4*c)-1)/2,z=(sqrt(1-4*c)+1)/2]  Show that z1+z2 = -1 (%i4) s:radcan(rhs(p+p)); (%o4) -1  Show that z2+z3=1 (%i6) s:radcan(rhs(p+p)); (%o6) 1  Reduced equation :  $z^{2}+z+c+1=0$ so standard coefficient names of single-variable quadratic polynomial $ax^{2}+bx+c$ are :  a = 1 b = 1 c = c+1  One can use the quadratic formula to find the exact solution ( roots of quadratic equation): $z_{1,2}={\frac {-b\pm {\sqrt {d}}}{2a}}={\frac {-1\pm {\sqrt {d}}}{2}}$ where : $d=b^{2}-4ac=1-4(c+1)=-4c-3$ Another formula can be used to solve square root of complex number " $s={\sqrt {d}}={\sqrt {r}}{\frac {d+r}{|d+r|}}$ where : $r=|d|$ so ###### Period 3 points In Maxima CAS kill(all); remvalue(z); f(z,c):=z*z+c; /* Complex quadratic map */ finverseplus(z,c):=sqrt(z-c); finverseminus(z,c):=-sqrt(z-c); /* */ fn(p, z, c) := if p=0 then z elseif p=1 then f(z,c) else f(fn(p-1, z, c),c); /*Standard polynomial F_p \, which roots are periodic z-points of period p and its divisors */ F(p, z, c) := fn(p, z, c) - z ; /* Function for computing reduced polynomial G_p\, which roots are periodic z-points of period p without its divisors*/ G[p,z,c]:= block( [f:divisors(p), t:1], /* t is temporary variable = product of Gn for (divisors of p) other than p */ f:delete(p,f), /* delete p from list of divisors */ if p=1 then return(F(p,z,c)), for i in f do t:t*G[i,z,c], g: F(p,z,c)/t, return(ratsimp(g)) )$

GiveRoots(g):=
block(
[cc:bfallroots(expand(%i*g)=0)],
cc:map(rhs,cc),/* remove string "c=" */
cc:map('float,cc),
return(cc)
)$compile(all); k:G[3,z,c]$
c:1;
s:GiveRoots(ev(k))$s;  ## Numerical methods ### Newton method #### find one root ##### algorithm For given ( and fixed = constant values): • period p • parameter c • initial aproximation $z_{0}$ = initial value of Newton iteration find one periodic point $z$ of complex quadratic polynomial $f$ :  $z:f^{p}(z)=z$ where  $f(z)=z^{2}+c$ Note that periodic point $z$ is only one of p points of periodic cycle. To find periodic point one have to solve equation ( = find a single zero of a function F): $f^{p}(z)-z=0$ Lets call left side of this equation (arbitrary name) function $F$ : $F_{p}(z)=f^{p}(z)-z$ Derivative of function $F$ with respect to variable z :  $F'_{p}(z)=f^{p'}(z)-1$ It can be checked using Maxima CAS: (%i1) f:z*z+c; (%o1) z^2 + c (%i2) F:f-z; (%o2) z^2 - z + c (%i3) diff(F,z,1); (%o3) 2 z - 1  Newton iteration is : $z_{n+1}=z_{n}-{\frac {F_{p}(z_{n})}{F'_{p}(z_{n})}}=z_{n}-{\frac {f_{p}(z_{n})-z_{n}}{f'_{p}(z_{n})-1}}=z_{n}-N_{p}(z_{n})$ gives sequence of $z_{n}$ points: $z_{0}$ $z_{1}=z_{0}-N_{p}(z_{0})$ $z_{2}=z_{1}-N_{p}(z_{1})$ $...$ $z_{n+1}=z_{n}-N_{p}(z_{n})$ where: • n is a number of Newton iteration • p is a period ( fixed positive integer) • $z_{0}$ is an initial aproximation of periodic point • $z_{n}$ is a final aproximation of periodic point • $f_{p}(z_{n})$ is computed using quadratic iteration with $z_{n}$ as a starting point • $f'_{p}(z_{n})$ is a value of first derivative . It computed using iteration • $N_{p}$ is called Newton function. Sometimes different and shorter notation is used : $N_{p}(z_{n})=z_{n}-{\frac {F_{p}(z_{n})}{F'_{p}(z_{n})}}$ ##### stoping criteria • $\epsilon _{succes}$ : predefined accuracy threshold ( succes = converged ) • $n_{max}$ : maximal number of iterations ( fail to converge) $\epsilon =|z_{n+1}-z_{n}|$ Usually one starts with :  $\epsilon _{succes}=10^{-16}$ for long double floating point number format  $n_{max}=10^{d}$ ##### precision Precision is proportional to the degree of polynomial F, which is proportional to the period p  $d=d(p)=2^{p-1}$ so  $period\to degree\to accuracy\to numberformat$ IEEE 754 - 2008 Common name C++ data type Base $b$ Precision $r$ Machine epsilon = $b^{-(r-1)}$ binary16 half precision short 2 11 2−10 ≈ 9.77e-04 binary32 single precision float 2 24 2−23 ≈ 1.19e-07 binary64 double precision double 2 53 2−52 ≈ 2.22e-16 extended precision _float80 2 64 2−63 ≈ 1.08e-19 binary128 quad(ruple) precision _float128 2 113 2−112 ≈ 1.93e-34 where : • precision = number of fraction's binary digits • accuracy < machine epsilon  precision = -log2(accuracy)  Relation between accuracy and degree d: • accuracy = 1/d • accuracy = 1/(2 d log d) Compare: ##### number of iterations ##### code ###### c /* gcc -std=c99 -Wall -Wextra -pedantic -O3 m.c -lm ./a.out */ #include <complex.h> #include <math.h> #include <stdio.h> #include <stdlib.h> const double pi = 3.141592653589793; static inline double cabs2(complex double z) { return (creal(z) * creal(z) + cimag(z) * cimag(z)); } /* newton function N(z) = z - (fp(z)-z)/fp'(z)) used for computing periodic point of complex quadratic polynomial f(z) = z*z +c */ complex double N( complex double c, complex double zn , int pMax, double er2){ complex double z = zn; complex double d = 1.0; /* d = first derivative with respect to z */ int p; for (p=0; p < pMax; p++){ d = 2*z*d; /* first derivative with respect to z */ z = z*z +c ; /* complex quadratic polynomial */ if (cabs(z) >er2) break; } z = zn - (z - zn)/(d - 1) ; return z; } /* compute periodic point of complex quadratic polynomial using Newton iteration = numerical method */ complex double GivePeriodic(complex double c, complex double z0, int period, double eps2, double er2){ complex double z = z0; complex double zPrev = z0; // prevoiuus value of z int n ; // iteration const int nMax = 64; for (n=0; n<nMax; n++) { z = N( c, z, period, er2); if (cabs2(z - zPrev)< eps2) break; zPrev = z; } return z; } int main() { const double ER2 = 2.0 * 2.0; // ER*ER const double EPS2 = 1e-100 * 1e-100; // EPS*EPS int period = 3; complex double c = -0.12565651+0.65720*I ; // https://commons.wikimedia.org/wiki/File:Fatou_componenets_3.png complex double z0 = 0.0 ; // initial aproximation for Newton method; usuallly crotical point complex double zp; // periodic point zp = GivePeriodic( c , z0, period, EPS2, ER2); // printf (" c = %f ; %f \n", creal(c), cimag(c)); printf (" period = %d \n", period); printf (" z0 = %f ; %f \n", creal(z0), cimag(z0)); printf (" zp = %.16f ; %.16f \n", creal(zp), cimag(zp)); return 0; }  ###### Maxima CAS /* batch file for Maxima CAS find periodic point of iterated function solve : c0: -0.965 + 0.085*%i$ //     period of critical orbit = 2
2						1 ( alfa ) repelling 						2						1 ( beta) repelling
[0.08997933216560844 %i - 0.972330689471866, 	0.0385332450804691 %i - 0.6029437025417168, 	(- 0.0899793321656081 %i) - 0.02766931052813359, 	1.602943702541716 - 0.03853324508046944 %i] // periodic z points
[1.952970301777068, 				1.208347498712429, 				0.1882750218384916, attr					3.206813573935709] // stability ( 1)
[0.3676955262170045, 				1.460103677644584, 				0.3676955262170032, attr					10.28365329797831] // stability (2)

Newton :
z0: 0 				 zn : - 0.08997933216560859 %i - 0.02766931052813337
z0: -0.86875-0.03125*%i 	 zn: 0.08997933216560858*%i-0.9723306894718666

(%o21) (-0.08997933216560859*%i)-0.02766931052813337
(%i22) zn2:Periodic(period,c0,(-0.86875)-0.03125*%i,ER,eps)
(%o22) 0.08997933216560858*%i-0.9723306894718666

*/

kill(all); /* if you want to delete all variables and functions. Careful: this does not completely reset the environment to its initial state. */
remvalue(z);
reset(); /*  Resets many (global) system variables */
ratprint : false; /* remove "rat :replaced  " */
display2d:false$/* display in one line */ mode_declare(ER, real)$
ER: 1e100$mode_declare(eps,real)$
eps:1e-100$mode_declare( c0, complex)$
c0: -0.965 + 0.085*%i$/* period of critical orbit = 2*/ mode_declare(z0, complex)$
z0:0$mode_declare(period, integer)$
period:2$/* ----------- functions ------------------------ */ /* https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/qpolynomials complex quadratic polynomial */ f(z,c):=z*z+c$

/* iterated map */

fn(p, z, c) :=
if p=0 then z
elseif p=1 then f(z,c)
else f(fn(p-1, z, c),c);

Fn(p,z,c):=  fn(p, z, c) - z$/* find periodic point using solve = symbolic method */ S(p,c):= ( s: solve(Fn(p,z,c)=0,z), s: map( rhs, s), s: float(rectform(s)))$

/* newton function : N(z) = z - (fp(z)-z)/(f'(z)-1) */

N(zn, c, iMax, _ER):=block(

[z, d], /* d = first derivative with respect to z */
/* mode_declare([z,zn, c,d], complex), */

d  : 1+0*%i,
z  : zn,

catch( for i:1 thru iMax step 1 do (

/*  print("i =",i,", z=",z, ", d = ",d),   */

d : float(rectform(2*z*d)), /* first derivative with respect to z */
z : float(rectform(z^2+c)), /* complex quadratic polynomial */

if ( cabs(z)> _ER)  /* z is escaping */
then  throw(z)

)), /* catch(for i */

if (d!=0) /* if not : divide by 0 error */
then  z:float(rectform( zn - (z - zn)/(d-1))),

return(z) /*  */

)$/* compute periodic point of complex quadratic polynomial using Newton iteration = numerical method */ Periodic(p, c, z0, _ER, _eps) :=block ( [z, n, temp], /* local variables */ if (p<1) then print("to small p ") else ( n:0, z: z0, /* */ temp:0, /* print("n = ", n, "zn = ", z ), */ catch( while (1) do ( z : N(z, c, p, _ER), n:n+1, /* print("n = ", n, "zn = ", z ), */ if(cabs(temp-z)< _eps) then throw(z), temp:z, if(n>64) then throw(z) ))), /* catch(while) */ return(z) )$

compile(all)$/* ------------------ run ---------------------------------------- */ zn1: Periodic(period, c0, z0, ER, eps); zn2 : Periodic(period, c0,-0.868750000000000-0.031250000000000*%i, ER, eps);  output: (%i21) zn1:Periodic(period,c0,z0,ER,eps) (%o21) (-0.08997933216560859*%i)-0.02766931052813337 (%i22) zn2:Periodic(period,c0,(-0.86875)-0.03125*%i,ER,eps) (%o22) 0.08997933216560858*%i-0.9723306894718666  #### find all roots Methods: • The iterated refinement Newton method • the Ehrlich-Aberth iteration ( Mpsolve ) Number: • number of all roots = degree of polynomial ( Fundamental Theorem of Algebra ) • degree = $2^{p}$ where p is a period period degree number of all roots 1 2^1 = 2 2 2 2^2=4 4 p d = 2^p d ##### test "Viete tests on the roots found: • if p is any polynomial of degree d and normalized (i.e. the leading coefficient is 1), then the sum of all roots of p must be equal to the negative of the coefficient of degree d − 1. • Moreover, the product of the roots must be equal to (−1)^d times the constant coefficient. If one root vanishes, then the product of the remaining d − 1 roots is equal to (−1)^(d−1) times the linear term. " $z_{1}+z_{2}+\dots +z_{d-1}+z_{d}=-{\dfrac {a_{d-1}}{a_{d}}}$ $z_{1}z_{2}\dots z_{d}=(-1)^{d}{\dfrac {a_{0}}{a_{d}}}$ "The sum of all roots should be 0 = the coefficient a of the degree d − 1 term of our polynomial. " $a_{d-1}=0$ $a_{d}=1$ The constant term is a0. Its value depends on the period (degree). ###### code /* batch file for Maxima CAS find all periodic point of iterated function using solve */ kill(all); /* if you want to delete all variables and functions. Careful: this does not completely reset the environment to its initial state. */ remvalue(all); reset(); /* Resets many (global) system variables */ ratprint : false; /* remove "rat :replaced " */ display2d:false$   /* display in one line */

mode_declare(c0, complex)$c0: -0.965 + 0.085*%i$ /*     period of critical orbit = 2*/
mode_declare(z, complex)$mode_declare(period, integer)$
period:2  $/* ----------- functions ------------------------ */ /* https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/qpolynomials complex quadratic polynomial */ f(z,c):=z*z+c$

/* iterated map */

fn(p, z, c) :=
if p=0 then z
elseif p=1 then f(z,c)
else f(fn(p-1, z, c),c);

Fn(p,z,c):=  fn(p, z, c) - z$/* find periodic point using solve = symbolic method */ S(p,c):= ( s: solve(Fn(p,z,c)=0,z), s: map( rhs, s), s: float(rectform(s)))$

/* gives a list of coefficients */
GiveCoefficients(P,x):=block
( [a, degree],
P:expand(P), /* if not expanded */
degree:hipow(P,x),
a:makelist(coeff(P,x,degree-i),i,0,degree));

GiveConstCoefficient(period, z, c):=(
[l:[]],
l :  GiveCoefficients( Fn(period,z,c), z),
l[length(l)]
)$/* https://en.wikipedia.org/wiki/Vieta%27s_formulas sum of all roots , here it should be zero */ VietaSum(l):=( abs(sum(l[i],i, 1, length(l))) )$

/*

it should be = (-1)^d* a0 = constant coefficient of Fn

*/

VietaProduct(l):=(
rectform(product(l[i],i, 1, length(l)))
)$compile(all); s: S(period, c0)$

vs: VietaSum(s);

vp : VietaProduct(s);

a0: GiveConstCoefficient( period, z,c0);


###### example results

The equation for period p=2 cycle is :

 z^4+2*c*z^2-z+c^2+c = 0


so d-1 = 3 coefficient is 0

$a_{d-1}=a_{3}=0$ $a_{0}=c^{2}+c$ Degree d :

$d=2^{p-1}=2^{2-1}=2^{2}=4$ The four roots for c= -0.965 + 0.085*i are :

[0.08997933216560844*%i-0.972330689471866, 0.0385332450804691*%i-0.6029437025417168, (-0.0899793321656081*%i)-0.02766931052813359, 1.602943702541716-0.03853324508046944*%i]


their sum is :

 $sum=6.938893903907228e-18\approx 0$ their product :

  $product=(-0.07904999999999943*\%i)-0.04100000000000011\approx a_{0}=c^{2}+c\approx (-0.07905*\%i)-0.04100000000000004$ So results have passed Viete's tests

# How to compute stability of periodic point ?

## stability of periodic points (orbit) - multiplier Critical orbit of discrete dynamical system based on complex quadratic polynomial. It tends to weakly attracting fixed point with abs(multiplier)=0.99993612384259

The multiplier (or eigenvalue, derivative) $m(f^{p},z_{0})=\lambda \,$ of a rational map $f\,$ iterated $p$ times at cyclic point $z_{0}\,$ is defined as:

$m(f^{p},z_{0})=\lambda ={\begin{cases}f^{p\prime }(z_{0}),&{\mbox{if }}z_{0}\neq \infty \\{\frac {1}{f^{p\prime }(z_{0})}},&{\mbox{if }}z_{0}=\infty \end{cases}}$ where $f^{p\prime }(z_{0})$ is the first derivative of$\ f^{p}$ with respect to $z\,$ at $z_{0}$ .

Because the multiplier is the same at all periodic points on a given orbit, it is called a multiplier of the periodic orbit.

The multiplier is:

• a complex number;
• invariant under conjugation of any rational map at its fixed point;
• used to check stability of periodic (also fixed) points with stability index $abs(\lambda ).\,$ A periodic point is

• attracting when $abs(\lambda )<1;\,$ • super-attracting when $abs(\lambda )=0;\,$ • attracting but not super-attracting when $0 • indifferent when $abs(\lambda )=1;\,$ • rationally indifferent or parabolic if $\lambda \,$ is a root of unity;
• irrationally indifferent if $abs(\lambda )=1\,$ but multiplier is not a root of unity;
• repelling when $abs(\lambda )>1.\,$ Periodic points

• that are attracting are always in the Fatou set
• that are repelling are in the Julia set;
• that are indifferent fixed points may be in one or the other. A parabolic periodic point is in the Julia set.

## code

### c

/*
gcc -std=c99 -Wall -Wextra -pedantic -O3 m.c -lm

c = -0.965000 ; 0.085000
period = 2
z0 = 0.000000 ; 0.000000
zp = -0.027669 ; -0.089979
m = 0.140000 ; 0.340000
Internal Angle = 0.187833

-----------------------------------------
for (int p = 1; p < pMax; p++){

if( cabs(c)<=2){
w =  MultiplierMap(c,p);
if ( iRadius <= 1.0) { // attractring
iAngle = GiveTurn(w);
period=p;
break;}
}   }

*/

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

const double pi = 3.141592653589793;

static inline double cabs2(complex double z) {
return (creal(z) * creal(z) + cimag(z) * cimag(z));
}

/* newton function : N(z) = z - (fp(z)-z)/f'(z)) */

complex double N( complex double c, complex double zn , int pMax, double er2){

complex double z = -zn;
complex double d = -1.0; /* d = first derivative with respect to z */
int p;

for (p=0; p < pMax; p++){

d = 2*z*d; /* first derivative with respect to z */
z = z*z +c ; /* complex quadratic polynomial */
if (cabs(z) >er2) break;

}

if ( cabs2(d) > 2)
z = zn - (z - zn)/d ;

return z;
}

/*
compute periodic point
using Newton iteration = numerical method

*/

complex double GivePeriodic(complex double c, complex double z0, int period, double eps2, double er2){

complex double z = z0;
complex double zPrev = z0; // prevoiuus value of z
int n ; // iteration
const int nMax = 64;

for (n=0; n<nMax; n++) {

z = N( c, z, period, er2);
if (cabs2(z - zPrev)< eps2) break;

zPrev = z; }

return z;
}

complex double AproximateMultiplierMap(complex double c, int period, double eps2, double er2){

complex double z;  // variable z
complex double zp ; // periodic point
complex double zcr = 0.0; // critical point
complex double d = 1;

int p;

// first find periodic point
zp =  GivePeriodic( c, zcr, period,  eps2, er2); // Find periodic point z0 such that Fp(z0,c)=z0 using Newton's method in one complex variable

// Find w by evaluating first derivative with respect to z of Fp at z0
if ( cabs2(zp)<er2) {

z = zp;
for (p=0; p < period; p++){
d = 2*z*d; /* first derivative with respect to z */
z = z*z +c ; /* complex quadratic polynomial */

}}
else d= 10000; //

return d;
}

complex double MultiplierMap(complex double c, int period, double eps2, double er2){

complex double w;
switch(period){
case 1: w = 1.0 - csqrt(1.0-4.0*c); 	break; // explicit
case 2: w = 4.0*c + 4; 			break; //explicit
default:w = AproximateMultiplierMap(c, period, eps2, er2); 	break; //

}

return w;

}

/*
gives argument of complex number in turns
https://en.wikipedia.org/wiki/Turn_(geometry)

*/

double GiveTurn( double complex z){
double t;

t =  carg(z);
t /= 2*pi; // now in turns
if (t<0.0) t += 1.0; // map from (-1/2,1/2] to [0, 1)
return (t);
}

int main() {

const double ER2  = 2.0 * 2.0; // ER*ER
const double EPS2 = 1e-100 * 1e-100; // EPS*EPS

int period = 3;
complex double c = -0.12565651+0.65720*I ; // https://commons.wikimedia.org/wiki/File:Fatou_componenets_3.png
complex double z0 = 0.0 ; // initial aproximation for Newton method; usuallly crotical point

complex double zp; // periodic point
complex double m; // multiplier

zp = GivePeriodic( c , z0, period, EPS2, ER2); //

m = MultiplierMap ( c, period,  EPS2, ER2 );

printf (" c = %f ; %f \n",  creal(c), cimag(c));
printf (" period = %d  \n",  period);
printf (" z0 = %f ; %f \n",  creal(z0), cimag(z0));
printf (" zp = %.16f ; %.16f \n",  creal(zp), cimag(zp));
printf (" m = %.16f ; %.16f \n",  creal(m), cimag(m));
printf (" Internal Angle = %.16f \n internal radius = %.16f \n", GiveTurn( m), cabs(m));

return 0;
}


### Maxima CAS

/*
Computes periodic points
z: z=f^n(z)
and its stability index ( abs(multiplier(z))

f(z)=z*z+c

batch file for Maxima cas
fraktal.republika.pl
*/

/* basic funtion  = monic and centered complex quadratic polynomial

*/
f(z,c):=z*z+c $/* iterated function */ fn(n, z, c) := if n=1 then f(z,c) else f(fn(n-1, z, c),c)$
/* roots of Fn are periodic point of  fn function */
Fn(n,z,c):=fn(n, z, c)-z $/* gives irreducible divisors of polynomial Fn[p,z,c] which roots are periodic points of period p */ G[p,z,c]:= block( [f:divisors(p),t:1], g, f:delete(p,f), if p=1 then return(Fn(p,z,c)), for i in f do t:t*G[i,z,c], g: Fn(p,z,c)/t, return(ratsimp(g)) )$

/* use :
roots:GiveRoots_bf(G[3,z,c]);

*/
GiveRoots_bf(g):=
block(

[cc:bfallroots(expand(g)=0)],
cc:map(rhs,cc),/* remove string "c=" */
return(cc)
)$/* -------------- main ----------------- */ load(cpoly); /* c should be inside Mandelbrot set to make attracting periodic points */ /* c:0.74486176661974*%i-0.12256116687665; center of period 3 component */ coeff:0.37496784+%i*0.21687214; /* -0.11+0.65569999*%i;*/ disp(concat("c=",string(coeff))); for p:1 step 1 thru 5 do /* period of z-cycle */ ( disp(concat("period=",string(p))), /* multiplier */ define(m(z),diff(fn(p,z,c),z,1)), /* G should be found when c is variable without value */ eq:G[p,z,c], /* here should be c as a symbol not a value */ c:coeff, /* put a value to a symbol here */ /* periodic points */ roots[p]:GiveRoots_bf(ev(eq)), /* ev puts value instead of c symbol */ /* multiplier of periodic points */ for z in roots[p] do disp(concat( "z=",string(float(z)),"; abs(multiplier(z))=",string(cabs(float(m(z)))) ) ), kill(c) )$

/*



results:

c:0.74486176661974*%i-0.12256116687665
periodic z-points:
supperattracting cycle (multiplier=0):
3.0879115871966895*10^-15*%i-1.0785250533182843*10^-14=0
0.74486176661974*%i-0.12256116687665=c
0.5622795120623*%i-0.66235897862236=c*c+c

=== repelling cycle
0.038593416246119-1.061217891258986*%i,
0.99358185708018-0.90887310643243*%i,

0.66294971900937*%i-1.247255127827274]
--------------------------------------------

for c:-0.11+0.65569999*%i;
abs(multiplier) = 0.94840556944351

z1a=0.17434140717278*%i-0.23080799291989,
z2a=0.57522120945524*%i-0.087122596659277,
z3a=0.55547045915754*%i-0.4332890929585,

abs(multiplier) =15.06988063154376
z1r=0.0085747730232722-1.05057855305735*%i,
z2r=0.95628667892607-0.89213756745704*%i,
z3r=0.63768304472883*%i-1.213641769411674