跳转到内容

分形/复平面的迭代/Mandelbrot 集/边界

来自 Wikibooks,开放世界中的开放书籍

Mandelbrot 集的边界包括:[1]

  • Mandelbrot 集的原始和卫星双曲分量的边界,包括点
    • 抛物线型(包括 1/4 和原始根,它们是具有有理外部角 = 可访问的 2 参数射线的着陆点)。
    • 西格尔型(具有无理外部角的唯一参数射线着陆)
    • 克雷默型(具有无理外部角的唯一参数射线着陆)
  • M 的边界,不包括具有点的双曲分量的边界
    • 不可重整化(米西乌列维奇型,具有有理外部角等)。
    • 有限重整化(米西乌列维奇型等)。
    • 无限重整化(费根鲍姆型等)。着陆在费根鲍姆点的外部射线的角度(以转数为单位)是无理数。
  • 非双曲分量的边界(我们认为它们不存在,但我们无法证明)。非双曲分量的边界也将是无限重整化的。


作为Mandelbrot 集 下的 单位圆的图像的边界


从 t 计算 c

[编辑 | 编辑源代码]
/*
gcc c.c -lm -Wall
./a.out
Root point between period 1 component and period 987 component  = c = 0.2500101310666710+0.0000000644946597
Internal angle (c) = 1/987
Internal radius (c) = 1.0000000000000000

*/

#include <stdio.h>
#include <math.h>		// M_PI; needs -lm also
#include <complex.h>






// numer of hyberolic components ( and it's centers ) of Mandelbrot set 
int GiveNumberOfCenters(int period){

	//int s = 0;
	int num = 1;
	
	int f; 
	int fMax = period-1; //sqrt(period); // https://stackoverflow.com/questions/11699324/algorithm-to-find-all-the-exact-divisors-of-a-given-integer
	
	
	
	if (period<1 ) {printf("input error: period should be positive integer \n"); return -1;}
	if (period==1) { return 1;}
		
	num = pow(2, period-1);
	
	// minus sum of number of center of all it's divisors (factors)
	for (f=1; f<= fMax; ++f){
	
		if (period % f==0)
			{num = num - GiveNumberOfCenters(f); }
	
	
	
	}
	
	
		
		
	


	return num ;

}


int ListNumberOfCenters(int period){

	int p=1;
	int pMax = period;
	int num=0;
	
	if (period ==1 || period==2) {printf (" for period %d there is only one component\n", period); return 0;}
	
	for (p=1; p<= pMax; ++p){
		num = GiveNumberOfCenters(p);
		printf (" for period %d there are %d components\n", p, num);
		}
		
	return 0;

}






/* 
   c functions using complex type numbers computes c from  component  of Mandelbrot set 
   input: 
   Period : positive integer
   n<d
   
   InternalRadius in [0.0, 1.0]
      
   */
complex double Give_c( int Period,  int n, int d , double InternalRadius )
{
  complex double c; 
  complex double w;  // point of reference plane  where image of the component is a unit disk 
 // alfa = ax +ay*i = (1-sqrt(d))/2 ; // result 
  double t; // InternalAngleInTurns
  
  if (Period<1 ) 
		{printf("input error: period should be positive integer \n"); return 10000.0;}
  
  
  
  t  = (double) n/d; 
  t = t * M_PI * 2.0; // from turns to radians
  
  w = InternalRadius*cexp(I*t); // map to the unit disk 
  
  switch ( Period ) // of component 
    {
    case 1: // main cardioid = only one period 1 component
      	c = w/2 - w*w/4; // https://wikibooks.cn/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set/boundary#Solving_system_of_equation_for_period_1
      	break;
    case 2: // only one period 2 component 
      	c = (w-4)/4 ; // https://wikibooks.cn/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set/boundary#Solving_system_of_equation_for_period_2
      	break;
      
    default: // period > 2 
      	printf("period = %d ; for higher periods : not works now, to do, use newton method \n", Period);
      	printf("for each Period >2 there are more then 1 components.\nHere are 2^(p-1) - s = %d components for period %d\n", GiveNumberOfCenters(Period), Period);
      	printf("to choose component use: it's center or external ray or angled internal address \n");
      	c = 10000.0; // bad value 
       
      break; }
      
      
  return c;
}



void PrintAndDescribe_c( int period,  int n, int d , double InternalRadius ){


	complex double c = Give_c(period, n, d, InternalRadius);
	
	printf("Root point between period %d component and period %d component  = c = %.16f%+.16f*I\t",period, d, creal(c), cimag(c));
	//printf("Internal radius (c) = %.16f\n",InternalRadius);
	printf("Internal angle (c) = %d/%d\n",n, d);
	



}

/*
https://stackoverflow.com/questions/19738919/gcd-function-for-c
The GCD function uses Euclid's Algorithm. 
It computes A mod B, then swaps A and B with an XOR swap.
*/

int gcd(int a, int b)
{
    int temp;
    while (b != 0)
    {
        temp = a % b;

        a = b;
        b = temp;
    }
    return a;
}


int ListAllBifurcationPoints(int period,  int d ){

	
	double InternalRadius = 1.0;
	// internal angle in turns as a ratio = n/d
	int n = 1;
	//int d = 987;
	
	if (period<1 ) 
		{printf("input error: period should be positive integer \n"); return 1;}
	if (period >2) 
		{printf("input error: not works now. TODO \n"); return 2;}
	// n/d = local angle in turns
 	for (n = 1; n < d; ++n ){
     		if (gcd(n,d)==1 )// irreducible fraction
       				{ PrintAndDescribe_c(period, n,d,InternalRadius); }
       		}       		
       	return 0;

}






int main (){

	int period = 1;
	double InternalRadius = 1.0;
	// internal angle in turns as a ratio = p/q
	int n = 1;
	int d = 987;
	complex double c;
	
	//c = Give_c(period, n,d , InternalRadius);
	//PrintAndDescribe_c(period, n,d , InternalRadius);
	
	ListAllBifurcationPoints(period,d);
	
	 ListNumberOfCenters(period);

	return 0;

}


从 c 计算 t

[编辑 | 编辑源代码]

描述参数 c 和内部角度 t 之间关系的函数

它用于计算主心形边界的 c 点

要从 c 计算 t,可以使用 Maxima CAS

(%i1) eq1:c = exp(%pi*t*%i)/2 -  exp(2*%pi*t*%i)/4;


                               %i %pi t     2 %i %pi t
                             %e           %e
(%o1)                    c = ---------- - ------------
                                 2             4
(%i2) solve(eq1,t);
             %i log(1 - sqrt(1 - 4 c))        %i log(sqrt(1 - 4 c) + 1)
(%o2) [t = - -------------------------, t = - -------------------------]
                    %pi                              %pi


所以

 f1(c):=float(cabs( -  %i* log(1 - sqrt(1 - 4* c))/%pi));
 f2(c):=float(cabs( -  %i* log(1 + sqrt(1 - 4* c))/%pi));



绘制边界

[编辑 | 编辑源代码]
收敛到 Mandelbrot 集边界的隐式多项式曲线 = 勒姆尼喀特

用于绘制 Mandelbrot 集边界的方法:[2]

如何绘制整个 M 集边界

[编辑 | 编辑源代码]

Jungreis 函数

[编辑 | 编辑源代码]
Mandelbrot 集的边界作为在 Jungreis 函数下单位圆的图像

描述

Python 代码

#!/usr/bin/env python
"""
    Python code by Matthias Meschede 2014
    http://pythology.blogspot.fr/2014/08/parametrized-mandelbrot-set-boundary-in.html
"""
import numpy as np
import matplotlib.pyplot as plt

nstore = 3000  #cachesize should be more or less as high as the coefficients
betaF_cachedata = np.zeros( (nstore,nstore))
betaF_cachemask = np.zeros( (nstore,nstore),dtype=bool)
def betaF(n,m):
    """
    This function was translated to python from
    http://fraktal.republika.pl/mset_jungreis.html
    It computes the Laurent series coefficients of the jungreis function
    that can then be used to map the unit circle to the Mandelbrot
    set boundary. The mapping of the unit circle can also
    be seen as a Fourier transform. 
    I added a very simple global caching array to speed it up
    """
    global betaF_cachedata,betaF_cachemask

    nnn=2**(n+1)-1
    if betaF_cachemask[n,m]:
        return betaF_cachedata[n,m]
    elif m==0:
        return 1.0
    elif ((n>0) and (m < nnn)):
        return 0.0
    else: 
        value = 0.
        for k in range(nnn,m-nnn+1):
            value += betaF(n,k)*betaF(n,m-k)
        value = (betaF(n+1,m) - value - betaF(0,m-nnn))/2.0 
        betaF_cachedata[n,m] = value
        betaF_cachemask[n,m] = True
        return value

def main():
    #compute coefficients (reduce ncoeffs to make it faster)
    ncoeffs= 2400
    coeffs = np.zeros( (ncoeffs) )
    for m in range(ncoeffs):
        if m%100==0: print '%d/%d'%(m,ncoeffs)
        coeffs[m] = betaF(0,m+1)

    #map the unit circle  (cos(nt),sin(nt)) to the boundary
    npoints = 10000
    points = np.linspace(0,2*np.pi,npoints)
    xs     = np.zeros(npoints)
    ys     = np.zeros(npoints)
    xs = np.cos(points)
    ys = -np.sin(points)
    for ic,coeff in enumerate(coeffs):
        xs += coeff*np.cos(ic*points)
        ys += coeff*np.sin(ic*points)
    
    #plot the function
    plt.figure()
    plt.plot(xs,ys)
    plt.show()

if __name__ == "__main__":
    main()

如何绘制双曲分量的边界

[编辑 | 编辑源代码]
周期为 1-6 的双曲分量的边界作为边界方程的解

绘制边界的几种方法

  • 求解边界方程
    • 使用 Brown、Stephenson、Jung 的方法。它仅适用于周期 7-8 [11]
    • 使用结果式 [12]
  • 使用 牛顿法 对边界进行参数化,该方法在分量中心附近有效 [13] [14]。该方法需要中心点,因此存在一些限制,
  • 边界扫描:通过检查每个像素,在检测到周期后检测边缘。这种方法很慢,但允许缩放。绘制“所有”分量
  • 针对给定 c 值的边界追踪。绘制单个分量。
  • Anne M. Burns 的假曼德布罗特集合:绘制主心形及其所有后代。不绘制迷你曼德布罗特集合。 [15]


"... 使用牛顿法绘制双曲分量的边界。也就是说,取您感兴趣的双曲分量中的一个点(那里有一个吸引循环),然后找到一个曲线,沿该曲线乘子的模数趋于 1。然后您将找到一个无差异参数。现在您可以同样改变乘子的幅角,同样使用牛顿法,并追踪这条曲线。在“尖点”附近需要格外小心。" Lasse Rempe-Gillen [16]


求解边界方程

[编辑 | 编辑源代码]
/* maxima CAS*/
kill(all);

f(z):=z^2+c;

d:diff(f(z),z,1);

r:1; /* only boundary */

/* system of equations */
e2:d = r*exp(2*%pi*%i*t);
e1:f(z) = z;

ss: solve([e1,e2],[z,c]);
s: ss[1]; 
s:s[2];
s:rhs(s);

/* change form */
x:realpart(s);
y:imagpart(s);

load(draw);
draw2d(
  line_type = solid,
  nticks= 500,
  parametric(x,y, t,0,2*%pi)
);


/* 
c functions using complex type numbers
computes c from  component  of Mandelbrot set */
complex double Give_c( int Period,  int p, int q , double InternalRadius )
{
  
  complex double w;  // point of reference plane  where image of the component is a unit disk 
  complex double c; // result 
  double t; // InternalAngleInTurns
  
  t  = (double) p/q; 
  t = t * M_PI * 2.0; // from turns to radians
  
  w = InternalRadius*cexp(I*t); // map to the unit disk 
  
  switch ( Period ) // of component 
    {
    case 1: // main cardioid = only one period 1 component
      c = w/2 - w*w/4; // https://wikibooks.cn/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set/boundary#Solving_system_of_equation_for_period_1
      break;
    case 2: // only one period 2 component 
      c = (w-4)/4 ; // https://wikibooks.cn/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set/boundary#Solving_system_of_equation_for_period_2
      break;
      // period > 2 
    default: 
    	printf("higher periods : to do, use newton method \n");
    	printf("for each q = Period of the Child component  there are 2^(q-1) roots \n");
      	c = 0.0; // 
       
      break; }
 return c;
}


定义周期为 n 的双曲分量边界的 2 个方程组

[编辑 | 编辑源代码]
  • 第一个定义周期点,
  • 第二个定义无差异轨道(周期点 的乘子等于 1)。

因为稳定指数 等于单位圆上的点的半径

所以可以将第二个方程改写为 [17]

它给出以下方程组


它可用于



以上 2 个方程组包含 3 个变量: 为常数,乘子 的函数)。必须删除一个变量才能求解。

边界是闭合曲线:心形或圆形。可以使用角度 (此处以圈数为单位,从 0 到 1)对边界上的点进行参数化。

进行评估后,可以将其代入上述系统,得到一个包含 2 个方程和 2 个变量 的系统。

现在可以求解。


对于周期

  • 1-3 可以使用符号方法求解,得到隐式(边界方程) 和显式函数(逆乘子映射) :
    • 1-2 很容易求解[18]
    • 3 可以使用“初等代数”(Stephenson)求解
  • >3 无法简化为显式函数,但
    • 可以简化为隐式解(边界方程) 并用数值方法求解
    • 可以使用 牛顿法求解非线性方程组 进行数值求解


示例

求解周期 1 的方程组
[edit | edit source]

以下是 Maxima 代码 



(%i4) p:1;
(%o4) 1
(%i5) e1:F(p,z,c)=z;
(%o5) z^2+c=z
(%i6) e2:m(p)=w;
(%o6) 2*z=w
(%i8) s:eliminate ([e1,e2], [z]);
(%o8) [w^2-2*w+4*c]
(%i12) s:solve([s[1]], [c]);
(%o12) [c=-(w^2-2*w)/4]
(%i13) define (m1(w),rhs(s[1]));
(%o13) m1(w):=-(w^2-2*w)/4

其中

  w:exp(2*%pi*%i*t);
  

(%i15) display2d:false;
(%o15) false
(%i16) 2*w;
(%o16) 2*%e^(2*%i*%pi*t)
(%i17) w*w;
(%o17) %e^(4*%i*%pi*t)

第二个方程只包含一个变量,可以消去该变量。由于边界方程很简单,因此很容易得到显式解

m1(w):=-(w^2-2*w)/4


所以

  m1(t) := block([w], w:exp(2*%pi*%i*t), return ((2*w-w^2)/4));
  g(t) := float(rectform(t));

数学方程

 
求解周期 2 的方程组
[edit | edit source]

以下是使用 Barton Willis 的 to_poly_solve 包的 Maxima 代码

(%i4) p:2;
(%o4) 2
(%i5) e1:F(p,z,c)=z;
(%o5) (z^2+c)^2+c=z
(%i6) e2:m(p)=w;
(%o6) 4*z*(z^2+c)=w
(%i7) e1:F(p,z,c)=z;
(%o7) (z^2+c)^2+c=z
(%i10) load(topoly_solver);
to_poly_solve([e1, e2], [z, c]);
(%o10) C:/PROGRA~1/MAXIMA~1.1/share/maxima/5.16.1/share/contrib/topoly_solver.mac
(%o11) [[z=sqrt(w)/2,c=-(w-2*sqrt(w))/4],[z=-sqrt(w)/2,c=-(w+2*sqrt(w))/4],[z=(sqrt(1-w)-1)/2,c=(w-4)/4],[z=-(sqrt(1-w)+1)/2,c=(w-4)/4]]
(%i12) s:to_poly_solve([e1, e2], [z, c]);
(%o12) [[z=sqrt(w)/2,c=-(w-2*sqrt(w))/4],[z=-sqrt(w)/2,c=-(w+2*sqrt(w))/4],[z=(sqrt(1-w)-1)/2,c=(w-4)/4],[z=-(sqrt(1-w)+1)/2,c=(w-4)/4]]
(%i14) rhs(s[4][2]);
(%o14) (w-4)/4
(%i16) define (m2 (w),rhs(s[4][2]));
(%o16) m2(w):=(w-4)/4

显式解 

m2(w):=(w-4)/4
求解周期 3 的方程组
[edit | edit source]

对于周期 3(以及更高周期),以前的方法没有得到结果(Maxima 代码) 

(%i14) p:3;
e1:z=F(p,z,c);
e2:m(p)=w;
load(topoly_solver);
to_poly_solve([e1, e2], [z, c]);
(%o14) 3
(%o15) z=((z^2+c)^2+c)^2+c
(%o16) 8*z*(z^2+c)*((z^2+c)^2+c)=w
(%i17) 
(%o17) C:/PROGRA~1/MAXIMA~1.1/share/maxima/5.16.1/share/contrib/topoly_solver.mac
`algsys' cannot solve - system too complicated.
#0: to_poly_solve(e=[z = ((z^2+c)^2+c)^2+c,8*z*(z^2+c)*((z^2+c)^2+c) = w],vars=[z,c])
-- an error.  To debug this try debugmode(true);

我使用 Robert P. Munafo[19] 编写的代码,该代码基于 Wolf Jung 的论文。

可以使用方程[20] 近似周期 3 分量

(%i1) z:x+y*%i;
(%o1)                              %i y + x
(%i2) w:asinh(z);
(%o2)                           asinh(%i y + x)
(%i3) realpart(w);
(%o3) 
                                                         2    2
          2    2     2      2  2 1/4     atan2(2 x y, - y  + x  + 1)      2
log((((- y  + x  + 1)  + 4 x  y )    sin(---------------------------) + y)
                                                      2
                                                        2    2
         2    2     2      2  2 1/4     atan2(2 x y, - y  + x  + 1)      2
 + (((- y  + x  + 1)  + 4 x  y )    cos(---------------------------) + x) )/2
                                                     2
(%i4) imagpart(w);
                                                                2    2
                 2    2     2      2  2 1/4     atan2(2 x y, - y  + x  + 1)
(%o4) atan2(((- y  + x  + 1)  + 4 x  y )    sin(---------------------------)
                                                             2
                                                              2    2
               2    2     2      2  2 1/4     atan2(2 x y, - y  + x  + 1)
     + y, ((- y  + x  + 1)  + 4 x  y )    cos(---------------------------) + x)
      

边界方程

[edit | edit source]

对上述系统求解的结果是 **边界方程**,


其中 是 **边界多项式**。

它定义了给定周期 的双曲分量的精确坐标。

这是一个隐式方程。

周期 1

可以轻松计算边界点 c

对于给定内角(旋转数)t,可以使用 Wolf Jung[21] 编写的代码计算周期 1 双曲分量(主心形)的边界点。

 t *= (2*PI); // from turns to radians
 cx = 0.5*cos(t) - 0.25*cos(2*t); 
 cy = 0.5*sin(t) - 0.25*sin(2*t);

周期 2

 t *= (2*PI);  // from turns to radians
 cx = 0.25*cos(t) - 1.0;
 cy = 0.25*sin(t);

求解边界方程

[edit | edit source]

对各种角度求解边界方程得到 边界点列表


可以使用 2 个复变量的牛顿法在特定分量中找到具有给定乘子值的点

使用方法类似

double _Complex z = 0;
double _Complex c = /* nucleus of component */;
int period = /* period of component */;
for (int t = 0; t < 360; ++t)
{
 m_d_interior(&z, &c, z, c, cexp(2 * M_PI * I * (t + 0.5) / 360), period, 16);
 /* plot c */
}

整个曼德勃罗集的面积

[编辑 | 编辑源代码]

分量的尺寸

[编辑 | 编辑源代码]


-- size-estimate.hs
-- http://www.fractalforums.com/mandelbrot-and-julia-set/some-questions-about-cardoids-circles-and-the-area/
import Data.Complex (Complex((:+)), polar)
import System.Environment (getArgs)

lambdabeta :: Int -> Complex Double -> (Complex Double, Complex Double)
lambdabeta period c = (lambda, beta)
  where
    zs = take period . iterate (\z -> z * z + c) $ 0
    lambdas = drop 1 . map (2 *) $ zs
    lambda = product lambdas
    beta = sum . map recip . scanl (*) 1 $ lambdas

sizeorient :: Int -> Complex Double -> (Double, Double)
sizeorient period c = polar . recip $ beta * lambda * lambda
  where
    (lambda, beta) = lambdabeta period c

main :: IO ()
main = do
  [p, x, y] <- getArgs
  let period = read p
      re     = read x
      im     = read y
      (size, orient) = sizeorient period (re :+ im)
  putStrLn $ "period = " ++ show period
  putStrLn $ "re     = " ++ show re
  putStrLn $ "im     = " ++ show im
  putStrLn $ "size   = " ++ show size
  putStrLn $ "orient = " ++ show orient

区分心形和伪圆

[编辑 | 编辑源代码]

方法是 区分心形和伪圆,描述在:A. Dolotin 著的通用曼德勃罗集。相关部分是第 5 节,特别是 5.8 方程。在该论文中, 在 5.1 方程中隐式定义,

然后 5.8 方程变为

当 epsilon 为

  • 接近 0 时,则为心形
  • 接近 1 时,则为圆形

代码

畸变和偏心率

[编辑 | 编辑源代码]
  • Md 的中心双曲分量的边界是具有 d - 1 个尖点的外摆线 [24][25]
  • 周期为 2 的分量的边界 [26]
  • 从 n 次分叉集的主分量中定位和计数卫星分量的分叉点 [27]
  • 一般 M-集周期区域中心的精确计算

参考文献

[编辑 | 编辑源代码]
  1. stackexchange:曼德勃罗集中的点分类
  2. mathoverflow:曼德勃罗集边界的参数化
  3. 边界扫描,作者:Robert P. Munafo,1993 年 2 月 3 日
  4. 杏仁面包主页
  5. 通过采样进行高效边界跟踪,作者:Alex Chen、Todd Wittman、Alexander Tartakovsky 和 Andrea Bertozzi
  6. http://www.robertnz.net/cx.htm 轮廓积分,作者:Robert Davies
  7. 维基百科中的 Jungreis 函数
  8. 使用 Jungreis 算法绘制 Mc
  9. 推特上 gro_tsen 关于曼德勃罗集的近似值,通过其边界的前 4096 个傅里叶系数:
  10. 图像和 Maxima CAS 源代码
  11. Mandelbrot 集成分量 : 带有描述和参考的图像
  12. Young Hee Geum,Kevin G. Hare : Groebner 基,结式和广义 Mandelbrot 集 - 使用因式分解来寻找不可约多项式,用于周期 2 及更高的周期
  13. 图像 : 使用牛顿法对边界进行参数化,靠近分量中心,并附带源代码
  14. Mark McClure "分岔集和临界曲线" - Mathematica 在教育和研究中的应用,第 11 卷,第 1 期 (2006)。
  15. Burns A M : 绘制逃逸:Mandelbrot 集中抛物线分岔的动画。数学杂志,第 75 卷,第 2 期 (2002 年 4 月),第 104-116 页
  16. mathoverflow : 如何有效地检测抛物线分量和西格尔圆盘分量,由 Lasse Rempe-Gillen 撰写
  17. 新闻组:gmane.comp.mathematics.maxima.general 主题:方程组 日期:2008-08-11 21:44:39 GMT
  18. Thayer Watkins : Mandelbrot 集的结构
  19. Brown 方法,由 Robert P. Munafo 撰写
  20. Mandelbrot 集的周期 3 双曲分量的参数化,由 Dante Giarrusso 和 Yuval Fisher 撰写
  21. Mandel:用于实数和复数动力学的软件,由 Wolf Jung 开发
  22. Mandelbrot 集和 Julia 集的组合学 - 1/n2 规则及其偏差
  23. fractalforums.org : 微型分形中的周期倍增
  24. 外摆线和布莱施克积
  25. Geum, Young Hee & Kim, Young. (2004). 度 n 分岔集主分量的外摆线边界。应用数学与计算杂志。16. 221-229. 10.1007/BF02936163.
  26. 度 3 分岔集周期 2 分量的参数边界,作者:Young Ik Kim,忠清数学学会杂志,第 16 卷,第 1 期,2003 年 6 月
  27. Geum, Young Hee & Kim, Young. (2006). 度 n 分岔集中卫星分量从主分量的分岔点定位和计数。应用数学与计算杂志。22. 339-350. 10.1007/BF02896483.
华夏公益教科书