Fractals/Iterations in the complex plane/periodic points

From Wikibooks, open books for an open world
Jump to navigation Jump to search

period[edit]

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

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 F[edit]

So equation for periodic points:

or in other notation

degree[edit]

Degree d of polynomial F :

How to find periodic points for given period?[edit]

Visual check[edit]

Visual check[2] of critical orbit helps to:

  • understand dynamics
  • find periodic points

Symbolic methods[edit]

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 divisors[edit]

  • math notation : [3]
  • 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 divisors[edit]

  • 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 points[edit]
Standard equations[edit]

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 equations[edit]

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
Computing periodic points[edit]

Definition in wikipedia [4]

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 :

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[edit]

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[5]

are :

 a = 1
 b = 1
 c = c+1

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


where :


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

where :

so

Period 3 points[edit]

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[edit]

Newton method[edit]

find one root[edit]

algorithm[edit]

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 criteria[edit]
  •  : predefined accuracy threshold ( succes = converged )[9]
  •  : maximal number of iterations ( fail to converge)

Usually one starts with :

 

for long double floating point number format

 
precision[edit]

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[10] 2 64 2−63 ≈ 1.08e-19
binary128 quad(ruple) precision _float128[10] 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:[11]

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

Compare:

number of iterations[edit]
code[edit]
c[edit]
/* 
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[edit]
/*

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[edit]

Methods:

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

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
p d = 2^p d


test[edit]

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

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

code[edit]
/*

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[edit]

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 reults have passed Viete's tests

How to compute stability of periodic point ?[edit]

stability of periodic points (orbit) - multiplier[edit]

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

A periodic point is[16]

  • attracting when
    • super-attracting when
    • attracting but not super-attracting when
  • indifferent when
    • rationally indifferent or parabolic 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.[17] A parabolic periodic point is in the Julia set.

code[edit]

c[edit]

/* 
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 CAS[edit]

/* 
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]

References[edit]

  1. [[::Periodic_point wikipedia : Periodic_point]]
  2. Period of critical orbit with slow dynamics
  3. Numerical Methods for Finding Periodic Orbits at Scholarpedia
  4. wikipedia : Periodic_points_of_complex_quadratic_mappings#Finite_fixed_points
  5. wikipedia : Quadratic function
  6. math.stackexchange question : how-to-find-the-solution-of-a-quadratic-equation-with-complex-coefficients
  7. wikipedia : Quadratic equation
  8. math.stackexchange question how-do-i-get-the-square-root-of-a-complex-number
  9. Newton's method in practice: finding all roots of polynomials of degree one million efficiently Dierk Schleicher, Robin Stoll
  10. a b Floating Types - Using the GNU Compiler Collection (GCC)
  11. 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
  12. wikipedia : MPSolve
  13. wikipedia : Vieta's formulas
  14. Newton's method in practice: finding all roots of polynomials of degree one million efficiently by Dierk Schleicher, Robin Stoll
  15. Alan F. Beardon, Iteration of Rational Functions, Springer 1991, ISBN 0-387-95151-2, p. 41
  16. Alan F. Beardon, Iteration of Rational Functions, Springer 1991, ISBN 0-387-95151-2, page 99
  17. Some Julia sets by Michael Becker