Fractals/Iterations in the complex plane/periodic points

How to describe a cycle ( periodic orbit)

  • period
  • stability of a periodic Orbit ( cycle)
  • (cyclic) permutation[1]
  • basin of attraction


periodEdit

Here period means period of iterated function (map).[2]

Point z has period p under f if :

or in other notation:


The smallest positive integer p satisfying the above is called the prime period or least period of the point z.

equation and polynomial FEdit

So equation for periodic points:

or in other notation

degreeEdit

Degree d of polynomial is the highest of the degrees of the polynomial's monomials (individual terms) with non-zero coefficients.

Degree of rational function is the higher value from the degree of numerator and denominator


In case of iterated quadratic polynomial  :

How to find periodic points for given period?Edit

 
Periodic points of f(z) = z*z-0.75 for period =6 as intersections of 2 implicit curves

Visual checkEdit

Visual check[3] of critical orbit helps to:

  • understand dynamics
  • find periodic points

Symbolic methodsEdit

Complex quadratic map

 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   which roots are periodic z-points of period p and its divisorsEdit

  • math notation :   [4]
  • ggplot[5]
  • Maxima CAS function :
F(p, z, c) := fn(p, z, c) - z ;

Function for computing reduced polynomial   which roots are periodic z-points of period p without its divisorsEdit

  • math definition :  
  • 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 pointsEdit
Standard equationsEdit

Equation for periodic point for period   is :

 

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  

These equations give periodic points for period p and its divisors. It is because these polynomials are product of lower period polynomials  :

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

 

 

 

 

 

Standard equations can be reduced to equations   giving periodic points for period p without its divisors.

See also : prime factors in wikipedia

Reduced equationsEdit

Reduced equation is

 

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



period degree Fp degree Gp
1 2^1 = 2 2
2 2^2=4 2
3 2^3 = 8 6
4 2^4 = 16 12
p d = 2^p
Computing periodic pointsEdit

Definition in wikipedia [6]

Fixed points ( period = 1 )Edit

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[1]));

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[1]+p[2]));
(%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 :


  
complex double GiveFixed(complex double c){
/* 

Equation defining fixed points : z^2-z+c = 0
	z*2+c = z
	z^2-z+c = 0

coefficients of standard form ax^2+ bx + c  
 	a = 1 , b = -1 , c = c
 
The discriminant d is 

	d=b^2- 4ac 
	d = 1 - 4c
	
 
 alfa =  (1-sqrt(d))/2 
*/

	complex double d = 1-4*c;
	complex double z = (1-csqrt(d))/2.0;
	return z;

}



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 pointsEdit

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[1]+p[2]));
(%o4) -1

Show that z2+z3=1

(%i6) s:radcan(rhs(p[3]+p[4]));
(%o6) 1

Reduced equation :

  

so standard coefficient names of single-variable quadratic polynomial[7]

 

are :

 a = 1
 b = 1
 c = c+1

One can use the quadratic formula[8] to find the exact solution ( roots of quadratic equation):[9]

 

where :

 

Another formula[10] can be used to solve square root of complex number "

 

where :

 

so

Period 3 pointsEdit

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 methodsEdit


Newton methodEdit

find one rootEdit

algorithmEdit

For given ( and fixed = constant values):

  • period p
  • parameter c
  • initial aproximation   = initial value of Newton iteration

find one periodic point   of complex quadratic polynomial  :

   

where

   

Note that periodic point   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):

 

Lets call left side of this equation (arbitrary name) function  :

 

Derivative of function   with respect to variable z :

  

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 :

 

gives sequence of   points:

 
 
 
 
 

where:

  • n is a number of Newton iteration
  • p is a period ( fixed positive integer)
  •   is an initial aproximation of periodic point
  •   is a final aproximation of periodic point
  •   is computed using quadratic iteration with   as a starting point
  •   is a value of first derivative . It computed using iteration
  •   is called Newton function. Sometimes different and shorter notation is used :  
stoping criteriaEdit
  •   : predefined accuracy threshold ( succes = converged )[11]
  •   : maximal number of iterations ( fail to converge)

 

Usually one starts with :

  

for long double floating point number format

  
precisionEdit

Precision is proportional to the degree of polynomial F, which is proportional to the period p

  

so

  
IEEE 754 - 2008 Common name C++ data type Base   Precision   Machine epsilon =  
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[12] 2 64 2−63 ≈ 1.08e-19
binary128 quad(ruple) precision _float128[12] 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:[13]

  • accuracy = 1/d
  • accuracy = 1/(2 d log d)

Compare:

number of iterationsEdit
codeEdit
cEdit
/* 
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 CASEdit
/*

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 rootsEdit

Methods:

  • The iterated refinement Newton method[13]
  • the Ehrlich-Aberth iteration ( Mpsolve[14] )

Number:

  • number of all roots = degree of polynomial ( Fundamental Theorem of Algebra )
  • degree =   where p is a period


period degree number of all roots
1 2^1 = 2 2
2 2^2=4 4
3 2^3 = 8 8
p d = 2^p d


testEdit

"Viete tests[15] 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. "

 

 

"The sum of all roots should be 0 = the coefficient a of the degree d − 1 term of our polynomial. "[16]

 

 

The constant term is a0. Its value depends on the period (degree).

codeEdit
/*

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 resultsEdit

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

 

 

Degree d :

 

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 :

  

their product :

   

So results have passed Viete's tests

How to compute stability of periodic point ?Edit

stability of periodic points (orbit) - multiplierEdit

 
Stability index of periodic points along horizontal axis
 
boundaries of regions of parameter plane with attracting orbit of periods 1-6
 
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)   of a rational map   iterated   times at cyclic point  is defined as:

 

where   is the first derivative of  with respect to   at  .

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;[17]
  • used to check stability of periodic (also fixed) points with stability index  

A periodic point is[18]

  • attracting when  
    • super-attracting when  
    • attracting but not super-attracting when  
  • indifferent ( neutral) when  
    • rationally indifferent or parabolic or weakly attractive[19] if   is a root of unity;
    • irrationally indifferent if   but multiplier is not a root of unity;
  • repelling when  

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.[20] A parabolic periodic point is in the Julia set.

codeEdit

cEdit

/* 
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 
 internal radius = 0.367696 
 
 
 -----------------------------------------
  for (int p = 1; p < pMax; p++){
      
      if( cabs(c)<=2){
      w =  MultiplierMap(c,p);
      iRadius = cabs(w);
      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 
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;
}

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 CASEdit

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

for complex quadratic polynomial :
f(z)=z*z+c

batch file for Maxima cas
Adam Majewski 20090604
fraktal.republika.pl
*/

/* basic funtion  = monic and centered complex quadratic polynomial 
http://en.wikipedia.org/wiki/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 :
load(cpoly); 
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

How to compute period of given periodic point ?Edit

How to find basin of attracting point ?Edit

What is the relation between Julia set and periodic point ?Edit

ReferencesEdit

  1. Very badly ordered cycles of interval maps by Sourav Bhattacharya, Alexander Blokh
  2. Periodic_point in wikipedia
  3. Period of critical orbit with slow dynamics
  4. Numerical Methods for Finding Periodic Orbits at Scholarpedia
  5. Fixed points in ggplot2 by Eric
  6. wikipedia : Periodic_points_of_complex_quadratic_mappings#Finite_fixed_points
  7. wikipedia : Quadratic function
  8. math.stackexchange question : how-to-find-the-solution-of-a-quadratic-equation-with-complex-coefficients
  9. wikipedia : Quadratic equation
  10. math.stackexchange question how-do-i-get-the-square-root-of-a-complex-number
  11. Newton's method in practice: finding all roots of polynomials of degree one million efficiently Dierk Schleicher, Robin Stoll
  12. a b Floating Types - Using the GNU Compiler Collection (GCC)
  13. a b Newton's method in practice II: The iterated refinement Newton method and near-optimal complexity for finding all roots of some polynomials of very large degrees by Robin Stoll, Dierk Schleicher
  14. wikipedia : MPSolve
  15. wikipedia : Vieta's formulas
  16. Newton's method in practice: finding all roots of polynomials of degree one million efficiently by Dierk Schleicher, Robin Stoll
  17. Alan F. Beardon, Iteration of Rational Functions, Springer 1991, ISBN 0-387-95151-2, p. 41
  18. Alan F. Beardon, Iteration of Rational Functions, Springer 1991, ISBN 0-387-95151-2, page 99
  19. Monireh Akbari, Maryam Rabii, Real cubic polynomials with a fixed point of multiplicity two, Indagationes Mathematicae, Volume 26, Issue 1, 2015, Pages 64-74, ISSN 0019-3577, https://doi.org/10.1016/j.indag.2014.06.001. (https://www.sciencedirect.com/science/article/pii/S0019357714000536) Abstract: The aim of this paper is to study the dynamics of the real cubic polynomials that have a fixed point of multiplicity two. Such polynomials are conjugate to fa(x)=ax2(x−1)+x, a≠0. We will show that when a>0 and x≠1, then |fan(x)| converges to 0 or ∞ and, if a<0 and a belongs to a special subset of the parameter space, then there is a closed invariant subset Λa of R on which fa is chaotic. Keywords: Multiplicity; Chaotic; Cubic polynomial; Schwarzian derivative
  20. Some Julia sets by Michael Becker