如何描述一个循环（周期轨道）

周期
周期轨道的稳定性
(循环) 置换^[1]
吸引域

周期

这里的周期是指迭代函数（映射）的周期。^[2]

如果点 z 在 f 下的周期为 p，则

z:\ f^{p}(z)=z

或用另一种表示法

z_{p}=z_{0}

满足上述条件的最小正整数 p 称为点 z 的素周期或最小周期。

方程和多项式 F

所以周期点的方程

f^{p}(z)=z

或用另一种表示法

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

度数

多项式的度数 d 是多项式中非零系数的单项式（单个项）的度数中最高的值。

有理函数的度数是分子和分母的度数中较高的那个值。

对于迭代的二次多项式

d=deg(f^{p})==2^{p}

如何为给定周期找到周期点？

一种用于计算一些多项式映射的精确 n 循环的算法

视觉检查

视觉检查^[3] 临界轨道有助于

理解动力学
找到周期点

周期为 1，动力学缓慢的轨道
周期为 3 的抛物型轨道
周期为 3 的轨道

符号方法

复二次映射

 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);

标准多项式Fp，其根是周期为p及其因子的周期z点

数学符号： $\ F_{p}(z,c)=f_{c}^{(p)}(z)-z$ ^[4]
ggplot^[5]
Maxima CAS 函数

F(p, z, c) := fn(p, z, c) - z ;

计算简化多项式Gp的函数，其根是周期为p的周期z点，但没有其因子的根

数学定义： $G_{p}(z,c)={\frac {F_{p}(z,c)}{\prod _{m|p,m<p}G_{m}(z,c)}}\,$
Maxima 函数

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))
)$

查找周期点的方程式

标准方程式

周期点的方程式对于周期 $p\,$ 是

$\ F_{p}(z,c)=0\,$

因此在 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

标准方程式的次数为 $2^{p}\,$

这些方程式给出了周期为p及其因子的周期点。这是因为这些多项式是较低周期多项式 $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)

所以

$\ 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}\,$

标准方程可以简化为方程 $G_{p}\,$ ，给出周期为p的周期点，但不包括其因子。

另请参见：维基百科中的素因子表

简化方程

简化方程为

$G_{p}(z,c)=0\,$

所以在 Maxima 中给出

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

周期	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

计算周期点

维基百科中的定义 ^[6]

不动点（周期 = 1）

定义不动点的方程

 z^2-z+c = 0

在 Maxima 中

定义 c 值

(%i1) c:%i;
(%o1) %i

计算不动点

(%i2) p:float(rectform(solve([z*z+c=z],[z])));
(%o2) [z=0.62481053384383*%i-0.30024259022012,z=1.30024259022012-0.62481053384383*%i]

找出哪些是排斥的

abs(2*rhs(p[1]));

证明不动点的总和为 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

绘制点

(%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]);

可以将点 z=1/2 添加到列表中

(%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]);

在 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;
}

周期为 2 的点

在 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]

证明 z1+z2 = -1

(%i4) s:radcan(rhs(p[1]+p[2]));
(%o4) -1

证明 z2+z3=1

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

简化方程

  $z^{2}+z+c+1=0$

所以单变量二次多项式的标准系数名称^[7]

ax^{2}+bx+c

是

 a = 1
 b = 1
 c = c+1

可以使用二次公式^[8] 找到精确解（二次方程的根）：^[9]

 $z_{1,2}={\frac {-b\pm {\sqrt {d}}}{2a}}={\frac {-1\pm {\sqrt {d}}}{2}}$

其中

 $d=b^{2}-4ac=1-4(c+1)=-4c-3$

另一个公式^[10] 可用于求解复数的平方根“

$s={\sqrt {d}}={\sqrt {r}}{\frac {d+r}{|d+r|}}$

其中

r=|d|

所以

周期为 3 的点

在 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;

数值方法

wolframalpha

牛顿法

找到一个根

算法

对于给定的（和固定 = 常数值）

周期 p
参数 c
初始近似值 $z_{0}$ = 牛顿迭代的初始值

寻找一个复二次多项式 $f$ 的周期点 $z$

   $z:f^{p}(z)=z$

其中

   $f(z)=z^{2}+c$

请注意，周期点 $z$ 只是周期循环中的 p 个点之一。

为了找到周期点，我们必须解方程（= 找到函数 F 的一个零点）

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

我们称此方程的左侧（任意名称）为函数 $F$

 $F_{p}(z)=f^{p}(z)-z$

函数 $F$ 关于变量 z 的导数

  $F'_{p}(z)=f^{p'}(z)-1$

可以使用 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

牛顿迭代是

 $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})$

给出点 $z_{n}$ 的序列

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})

其中

n 是牛顿迭代的次数
p 是周期（固定的正整数）
$z_{0}$ 是周期点的初始近似值
$z_{n}$ 是周期点的最终近似值
$f_{p}(z_{n})$ 使用以 $z_{n}$ 为起点的二次迭代计算
$f'_{p}(z_{n})$ 是第一次导数的值。使用迭代计算
$N_{p}$ 称为牛顿函数。有时使用不同的、更短的符号： $N_{p}(z_{n})=z_{n}-{\frac {F_{p}(z_{n})}{F'_{p}(z_{n})}}$

停止标准

$\epsilon _{succes}$ : 预定义的精度阈值 (succes = 收敛)^[11]
$n_{max}$ : 最大迭代次数 (无法收敛)

$\epsilon =|z_{n+1}-z_{n}|$

通常从

  $\epsilon _{succes}=10^{-16}$

对于长双精度浮点数格式

  $n_{max}=10^{d}$

精度

精度与多项式 F 的次数成正比，该次数与周期 p 成正比

  $d=d(p)=2^{p-1}$

所以

  $period\to degree\to accuracy\to numberformat$

IEEE 754 - 2008	常用名称	C++ 数据类型	基数 $b$	精度 $r$	机器精度 = $b^{-(r-1)}$
二进制16	半精度	short	2	11	2⁻¹⁰ ≈ 9.77e-04
二进制32	单精度	float	2	24	2⁻²³ ≈ 1.19e-07
二进制64	双精度	double	2	53	2⁻⁵² ≈ 2.22e-16
	扩展精度	_float80^[12]	2	64	2⁻⁶³ ≈ 1.08e-19
二进制128	四(重)精度	_float128^[12]	2	113	2⁻¹¹² ≈ 1.93e-34

其中

精度 = 小数部分的二进制位数
精度 < 机器精度

 precision =  -log2(accuracy)

精度和次数 d 之间的关系:^[13]

精度 = 1/d
精度 = 1/(2 d log d)

比较

使用 Maxima CAS 计算中心点的精度
缩放的精度

迭代次数

代码

来自 mandel 的 cpp

c

来自 mandelbrot-numerics 库的 m-attractor

/* 
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 critical 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://wikibooks.cn/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);

输出

(%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

找到所有根

方法

迭代细化牛顿法^[13]
Ehrlich-Aberth 迭代 (Mpsolve^[14])

数量

所有根的数量 = 多项式的次数 (代数基本定理)
次数 = $2^{p}$ 其中 p 是一个周期

周期	度数	所有根的数量
1	2^1 = 2	2
2	2^2=4	4
3	2^3 = 8	8
p	d = 2^p	d

测试

"维埃特公式^[15] 对找到的根进行测试

如果 p 是任何次数为 d 的多项式，并且是归一化的 (即最高次项系数为 1)，则 p 的所有根的和必须等于 d-1 次项系数的负值。
此外，根的乘积必须等于 (-1)^d 乘以常数项。如果一个根为零，则剩余 d-1 个根的乘积等于 (-1)^(d-1) 乘以一次项。"

$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}}}$

"所有根的总和应为 0 = 多项式的 d − 1 次项的系数 a。"^[16]

$a_{d-1}=0$

$a_{d}=1$

常数项为 a0。它的值取决于周期（次数）。

代码

/*

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://wikibooks.cn/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);

示例结果

周期 p=2 循环的方程是

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

所以 d-1 = 3 系数为 0

$a_{d-1}=a_{3}=0$

$a_{0}=c^{2}+c$

次数 d

 $d=2^{p-1}=2^{2-1}=2^{2}=4$

c= -0.965 + 0.085*i 的四个根是

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

它们的总和是

  $sum=6.938893903907228e-18\approx 0$

它们的乘积是

   $product=(-0.07904999999999943*\%i)-0.04100000000000011\approx a_{0}=c^{2}+c\approx (-0.07905*\%i)-0.04100000000000004$

所以结果通过了韦达定理的检验

如何计算周期点的稳定性？

周期点（轨道）的稳定性 - 乘子

乘子（或特征值，导数) $m(f^{p},z_{0})=\lambda \,$ 是有理映射 $f\,$ 在循环点 $z_{0}\,$ 处迭代 $p$ 次时的定义如下

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}}

其中 $f^{p\prime }(z_{0})$ 是 $\ f^{p}$ 关于 $z\,$ 在 $z_{0}$ 处的导数。

因为乘数在给定轨道上的所有周期点都是相同的，所以它被称为周期轨道的乘数。

乘数是

一个复数；
在任何有理映射在其不动点处的共轭下是不变的；^[17]
用于检查周期点（也称为不动点）的稳定性，稳定性指数为 $abs(\lambda ).\,$

一个周期点^[18]

当 $abs(\lambda )<1;\,$ $abs(\lambda )<1;\,$ 时是吸引的；
- 当 $abs(\lambda )=0;\,$ 时是超吸引的；
- 当 $0<abs(\lambda )<1;\,$ 时是吸引的，但不是超吸引的；
当 $abs(\lambda )=1;\,$ $abs(\lambda )=1;\,$ 时是中性的（无差别的）；
- 如果 $\lambda \,$ 是单位根，则它是理性无差别的或抛物线的或弱吸引的^[19]；
- 如果 $abs(\lambda )=1\,$ 但乘数不是单位根，则它是无理性无差别的；
当 $abs(\lambda )>1.\,$ 时是排斥的。

周期点

如果是吸引的，则始终位于法图集；
如果是排斥的，则位于朱利亚集；
如果是无差别的不动点，则可能位于两者之一。^[20] 抛物线周期点位于朱利亚集中。

代码

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 critical 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 function  = 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)
)$

/*

结果

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

如何计算给定周期点的周期？

如何计算真实周期？

Robert Munafo 在 Mu-Ency 中描述的多边形迭代法
Knighty 的泰勒球迭代法

两者都可以在一个区域中找到最低周期的核。如果该区域中没有核（要么 100% 外部，要么 100% 内部不围绕核），我不知道会发生什么。

多边形方法可以适应前周期 Misiurewicz 点。泰勒球方法可能也可以。球方法的一个版本可以适应使用雅可比矩阵的燃烧船等，而多边形方法有时会因早期折叠而失败。

int GivePeriod(complex double c){




	if (cabs(c)>2.0) {return 0;} // exterior
	if (cabs(1.0 - csqrt(1.0-4.0*c))<=1.0 ) {return 1;} // main cardioid
	if (cabs(4.0*c + 4)<=1.0){return 2;} // period 2 component
	
	int period =  GivePeriodByIteration(c); 
	
	if ( period < 0) // last chance
		{ 
			iUnknownPeriod +=1;
			//fprintf(stderr,"unknown period = %d < 0 c = %.16f %+.16f \n", period, creal(c), cimag(c) ); // print for analysis
			//period = m_d_box_period_do(c, 0.5, iterMax_LastIteration);
			//fprintf(stderr,"box period = %d  \n", period);
			
		}
	
	// period > 0 means is periodic 
	// period = 0 means is not periodic = exterior = escaping to infinity
	// period < 0 means period not found, maybe increase global variable iterMax_Period ( see local_setup)

	return period; 

}

marcm200

A heuristics to get to a point where the period detection might work correctly could also be: Compute the rolling derivative from the start on, i.e 2^n*product over all iterates. If the orbit approaches an attracting cycle of  length p, cyclic portions of that overall product will become smaller than one in magnitude. And repeatedly multiplying those values will eventually let the |rolling derivative| fall below any ever so small threshold. So if one  were to define beforehand such a small value like 10^-40, if the rolling derivative is first below, then doing a periodicity check by number equaltiy up to an epsilon should give the correct period.  I only use it with unicritical polynomials so chances are that during iteration one does not fall onto a critical point and then the rolling derivative gets stuck at 0 forever (so for z^2+c, starting at c is usually feasible, as +-sqrt(-c) is then almost always irrationally complex and not representable in software. However, with bad luck, its approximation might still lead to computer-0 within machine-epsilon).

( marcm200)

Knighty 的泰勒球迭代法

如 knightly 所建议的那样，找到 SA 多项式的根是一个完全不同的方法。您使用哪种多项式求根算法？

ball interval arithmetic" is also called "midpoint-radius interval".

Divide and conquer, using ball interval arithmetic to isolate the roots and Newton-Raphson to refine them.
My approach to Ball interval arithmetic is  very simple:
To a given number z we associate a disc representing an estimation of the errors. A ball interval may be represented like this: [z] = z + r [E]

Here z is the center of the ball, r is the radius and [E] is the error symbol representing the unit disc centred at 0.

Being not a mathematician, be aware that these definitions are very very sloppy. In particular, I don't care about rounding stuff.

One can also think about [E] as a function in some "mute" variables. for example, for complex numbers we can write:
[E] = f(a,b) = a * exp(i *b), where a is in [0,1] and b any real number.
it can also be written like this:
[E]= [0,1] * exp(i * [-inf, +inf])

A nice property of [E] is that raising it to a (positive) power will yield [E].
Another one is: [E] = -[E]
... Another one: z[E] = |z| [E]

Some basic operations on ball interval are:
-Addition:
  [v] + [w] = v + rv [E] + w + rw [E] = (v+w) + (rv+rw) [E]
-Subtraction:
  [v] - [w] = (v-w) + (rv+rw) [E]; because [E] will absorb the '-' sign.
-Multiplication:
 [v] [w] = (v + rv [E]) (w + rw [E]) = v*w + v*rw [E] + w*rv [E] + rv*rw [E]? = v*w + (|v|*rw + |w|*rv + rv*rw) [E]; because [E]? = [E]
The division is more complicated: It is not always defined (when 0 is inside the ball) and needs special care.
-Raising to a power:
use successive multiplications... but that may give too large balls... there are other methods.

With these operations one can evaluate polynomials using Horner method for example. It can also be used to iterate zn+1 = zn? + c for example and be used to find nucleus and its period instead of the box method.

您可以查看以下程序的代码：

查看此笔记

如何找到吸引点的盆地？

牛顿-拉夫森分形

朱利亚集和周期点之间有什么关系？

另请参阅

wikibooks 中的 critical_orbit

代码

关于使用自相关和函数逼近检测时间序列数据中的周期和重复模式的一种方法作者：Tim Breitenbach a b，Bartosz Wilkusz b，Lauritz Rasbach b，Patrick Jahnke
- Lauritz R：python 中的周期检测
- python 周期检测

参考文献

[1] Sourav Bhattacharya 和 Alexander Blokh 撰写的《区间映射的极度无序循环》

[2] 维基百科中的周期点

[3] 缓慢动力学下的临界轨道周期

[4] Scholarpedia 中查找周期轨道的数值方法

[5] Eric 撰写的 ggplot2 中的固定点

[6] 维基百科：复二次映射的周期点#有限固定点

[7] 维基百科：二次函数

[8] th.stackexchange 问题：如何找到具有复系数的二次方程的解

[9] 维基百科：二次方程

[10] th.stackexchange 问题：如何求复数的平方根

[11] Dierk Schleicher 和 Robin Stoll 撰写的《牛顿法在实践中：有效地查找一百万次方多项式的所有根》

[xfloat-12] 浮点类型 - 使用 GNU 编译器集合 (GCC)

[arxiv.org-13] Robin Stoll 和 Dierk Schleicher 撰写的《牛顿法在实践中 II：迭代细化牛顿法和查找一些非常大次数的多项式的所有根的近乎最优复杂度》

[14] 维基百科：MPSolve

[15] 维基百科：韦达定理

[16] Dierk Schleicher 和 Robin Stoll 撰写的《牛顿法在实践中：有效地查找一百万次方多项式的所有根》

[17] Alan F. Beardon 撰写的《有理函数的迭代》，Springer 1991，ISBN 0-387-95151-2，第 41 页

[18] Alan F. Beardon 撰写的《有理函数的迭代》，Springer 1991，ISBN 0-387-95151-2，第 99 页

[19] Monireh Akbari 和 Maryam Rabii 撰写的《具有双重固定点的实三次多项式》，《Indagationes Mathematicae》，第 26 卷，第 1 期，2015 年，第 64-74 页，ISSN 0019-3577，https://doi.org/10.1016/j.indag.2014.06.001. (https://www.sciencedirect.com/science/article/pii/S0019357714000536) 摘要：本文旨在研究具有双重固定点的实三次多项式的动力学。此类多项式与 fa(x)=ax2(x−1)+x，a≠0 共轭。我们将证明当 a>0 且 x≠1 时，|fan(x)| 收敛于 0 或 ∞，如果 a<0 且 a 属于参数空间的某个特殊子集，则存在 R 上的一个封闭不变子集 Λa，在该子集上 fa 是混沌的。关键词：重数；混沌；三次多项式；施瓦茨导数

[20] Michael Becker 撰写的《一些 Julia 集》

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

周期

方程和多项式 F

度数

如何为给定周期找到周期点？

视觉检查

符号方法

标准多项式Fp，其根是周期为p及其因子的周期z点

计算简化多项式Gp的函数，其根是周期为p的周期z点，但没有其因子的根

查找周期点的方程式

标准方程式

简化方程

计算周期点

不动点（周期 = 1）

周期为 2 的点

周期为 3 的点

数值方法

牛顿法

找到一个根

算法

停止标准

精度

迭代次数

代码

c

Maxima CAS

找到所有根

测试

代码

示例结果

如何计算周期点的稳定性 ？

周期点（轨道）的稳定性 - 乘子

代码

c

Maxima CAS

如何计算给定周期点的周期？

marcm200

Knighty 的泰勒球迭代法

如何找到吸引点的盆地？

朱利亚集和周期点之间有什么关系？

另请参阅

代码

参考文献

如何计算周期点的稳定性？