Fractals/Iterations in the complex plane/periodic points
period
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
So equation for periodic points:
or in other notation
degree
Degree d of polynomial F :
How to find periodic points for given period?
Visual check
Visual check[2] of critical orbit helps to:
- understand dynamics
- find periodic points
-
period 1 orbit with slow dynamics
-
period 3 parabolic orbit
-
period 3 orbit
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 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 which roots are periodic z-points of period p without its divisors
- 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
Standard equations
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
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 points
Definition in wikipedia [5]
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[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 :
<source lang=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[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[6]
are :
a = 1 b = 1 c = c+1
One can use the quadratic formula[7] to find the exact solution ( roots of quadratic equation):[8]
where :
Another formula[9] can be used to solve square root of complex number "
where :
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 = 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
- : predefined accuracy threshold ( succes = converged )[10]
- : maximal number of iterations ( fail to converge)
Usually one starts with :
for long double floating point number format
precision
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[11] | 2 | 64 | 2−63 ≈ 1.08e-19 | |
binary128 | quad(ruple) precision | _float128[11] | 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:[12]
- accuracy = 1/d
- accuracy = 1/(2 d log d)
Compare:
- precision for computing centers with Maxima CAS
- precision for zoom
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:
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 |
test
"Viete tests[14] 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. "[15]
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
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 ?
stability of periodic points (orbit) - multiplier
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;[16]
- used to check stability of periodic (also fixed) points with stability index
A periodic point is[17]
- 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.[18] 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
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
/* 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 ?
How to find basin of attracting point ?
What is the relation between Julia set and periodic point ?
- Parameter-Plane Dynamics of Fixed Points and Their Preimages For Standard Quadratic Julia Sets
- Quadratic Julia sets and periodic cycles
References
- ↑ Periodic_point in wikipedia
- ↑ Period of critical orbit with slow dynamics
- ↑ Numerical Methods for Finding Periodic Orbits at Scholarpedia
- ↑ Fixed points in ggplot2 by Eric
- ↑ wikipedia : Periodic_points_of_complex_quadratic_mappings#Finite_fixed_points
- ↑ wikipedia : Quadratic function
- ↑ math.stackexchange question : how-to-find-the-solution-of-a-quadratic-equation-with-complex-coefficients
- ↑ wikipedia : Quadratic equation
- ↑ math.stackexchange question how-do-i-get-the-square-root-of-a-complex-number
- ↑ Newton's method in practice: finding all roots of polynomials of degree one million efficiently Dierk Schleicher, Robin Stoll
- ↑ a b Floating Types - Using the GNU Compiler Collection (GCC)
- ↑ 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
- ↑ wikipedia : MPSolve
- ↑ wikipedia : Vieta's formulas
- ↑ Newton's method in practice: finding all roots of polynomials of degree one million efficiently by Dierk Schleicher, Robin Stoll
- ↑ Alan F. Beardon, Iteration of Rational Functions, Springer 1991, ISBN 0-387-95151-2, p. 41
- ↑ Alan F. Beardon, Iteration of Rational Functions, Springer 1991, ISBN 0-387-95151-2, page 99
- ↑ Some Julia sets by Michael Becker