跳轉到內容

分形/複平面的迭代/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://wikibook.tw/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://wikibook.tw/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://wikibook.tw/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://wikibook.tw/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.
華夏公益教科書