跳至內容

分形/複平面迭代/曼德勃羅集/中心

來自華夏公益教科書,開放書籍,開放世界

名稱 (同義詞)

[編輯 | 編輯原始碼]

中心 是一個雙曲分量 H 的引數 (或引數平面的點),使得對應的週期軌道具有乘子= 0。"[2]

初始域

[編輯 | 編輯原始碼]
  • "所有 M 的任何週期的中心都包含在 M 的內部,因此圓 C = {z ∈ C: |z − 0.75| = 2} 包圍著所有中心。"
  • "可以利用 M 的實對稱性,只使用具有非負虛部的起點;"[3]

中心的週期 = 臨界軌道的週期 = 雙曲分量的週期


朱利亞集

[編輯 | 編輯原始碼]

它是繪製朱利亞集最簡單的情況。

動力學平面由兩個超吸引盆組成

  • 外部 (無窮大是一個超吸引不動點)
  • 內部 (z=0 是超吸引週期點之一)


參見通用類別

曼德勃羅集元件的數量

[編輯 | 編輯原始碼]

方程

其中

  • N(p) 是給定週期 p 的曼德勃羅集的所有元件的數量。
  • p 是週期
  • s 是其各個分量在小於週期 p 的週期每個因數的元件數量之和[4]

所有值都是正整數。它是序列 A000740

Maxima CAS 中的程式

N(p) := block( 
[a:2^(p-1),f:0], 
for f in divisors(p) do 
   if f < p then a : a - N(f), 
return(a)
);

您可以執行它

for p:1 thru 50 step 1 do display(N(p));
N(1)=1
N(2)=1
N(3)=3
N(4)=6
N(5)=15
N(6)=27
N(7)=63
N(8)=120
N(9)=252
N(10)=495
N(11)=1023
N(12)=2010
N(13)=4095
N(14)=8127
N(15)=16365
N(16)=32640
N(17)=65535
N(18)=130788
N(19)=262143
N(20)=523770
N(21)=1048509
N(22)=2096127
N(23)=4194303
N(24)=8386440
N(25)=16777200
N(26)=33550335
N(27)=67108608
N(28)=134209530
N(29)=268435455
N(30)=536854005
N(31)=1073741823
N(32)=2147450880
N(33)=4294966269
N(34)=8589869055
N(35)=17179869105
N(36)=34359605280
N(37)=68719476735
N(38)=137438691327
N(39)=274877902845
N(40)=549755289480
N(41)=1099511627775
N(42)=2199022198821
N(43)=4398046511103
N(44)=8796090925050
N(45)=17592186027780
N(46)=35184367894527
N(47)=70368744177663
N(48)=140737479934080
N(49)=281474976710592
N(50)=562949936643600

sum(N(p),p,1,50);

結果

 1125899873221781
 
 

還有 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;

}


int main (){

	int period = 15;
	
	
	
	 ListNumberOfCenters(period);

	return 0;

}


另請參見

如何計算中心?

[編輯 | 編輯原始碼]

方法

  • 給定週期 的所有中心
    • 求解 [5]
      • 以圖形形式顯示結果
      • 使用數值方法[6][7]
        • “快速而骯髒”的演算法:檢查是否,則用顏色 n 給 c 點著色。這裡 n 是吸引軌道週期,eps 是吸引點周圍圓的半徑 = 數值計算的精度。
        • 找到多項式的所有根
          • 迭代細化牛頓法[8]
          • Ehrlich-Aberth 迭代(Mpsolve[9]
        • 原子域
        • Myrberg 方法[10][11]
      • 區間算術和分治策略[12]
    • “搜尋是使用網格的四叉樹細分進行的,這會在根密度更高的區域提高搜尋解析度。在一個搜尋區域內,首先以中等解析度繪製函式幅度的圖,然後使用 640 位定點數學細化區域性最小值點。由於乘以非常小的數字,根搜尋需要目標結果的兩倍精度。隨著根的新增,四叉樹進一步細分,搜尋繼續。由於對稱性,只搜尋上半平面。

在更高的週期,四叉樹深度被有意地限制,以減少執行時間,因為根的數量變得難以管理。” James Artherton

  • 給定週期,給定點 附近的一箇中心
    • 求解牛頓法 的一些步驟

定義中心的方程組

[edit | edit source]
週期為 10 的雙曲分量的中心
使用 不可約多項式 計算的週期為 1-10 的雙曲分量的中心

雙曲分量的中心是一個引數平面的點,對於該點,週期性的 z 迴圈是超吸引的。它給出了一個定義週期為 的雙曲分量的中心的兩個方程組

  • 首先定義週期點
  • 其次使 點超吸引。

有關詳細資訊,請參見 定義

“這些多項式具有整數係數。它們可以透過關於度的遞推關係獲得。令 Pd 為度數為 d 的多項式。我們有” [13]

 
 

求解方程組

[edit | edit source]

當臨界點 在此迴圈中時,可以求解第二個方程

為了解決系統,將 代入第一個方程。

定義中心的方程
[編輯 | 編輯原始碼]

得到方程:

分量的中心是上述方程的根。

因為 ,我們可以從這些方程中去掉 z

  • 對於週期 1:z^2+c=z 和 z=0 ,因此
  • 對於週期 2: (z^2+c)^2 +c =z 和 z=0 ,因此
  • 對於週期 3: ((z^2+c)^2 +c)^2 +c = z 和 z=0 ,因此

以下是計算上述函式的 Maxima 函式:

P[n]:=if n=1 then c else P[n-1]^2+c;
簡化後的定義中心方程
[編輯 | 編輯原始碼]

上述方程的根集包含週期為 p 及其因數的中心的集合。這是因為

其中

  • = 格利森週期為 m 的多項式 [14] = 的不可約因式 [15] [16] [1]
  • 使用迭代乘法的大寫 Pi 符號
  • 在這裡的意思是:對於所有 因數(從 1 到 p)。參見因數表

例如:

因此,可以使用以下方法找到不可約多項式

這是用於計算 的 Maxima 函式

GiveG[p]:=
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(Fn(p,0,c)),
for i in f do t:t*GiveG[i],
g: Fn(p,0,c)/t,
return(ratsimp(g))
)$

以下是週期為 1 到 10 的 的階數表,以及計算這些函式的係數所需的精度(只要使用未展開的形式 進行運算,根就可以用更低的精度來計算)。

週期
1 1 1 16
2 2 1 16
3 4 3 16
4 8 6 16
5 16 15 16
6 32 27 16
7 64 63 32
8 128 120 64
9 256 252 128
10 512 495 300
11 1024 1023 600

這裡

  • fpprec 是 bigfloat 數值的算術運算中有效的小數位數[17]

以下是最大系數表。

週期
1 0 1
2 0 1
3 1 2
4 1 3
5 3 116
6 4 5892
7 11 17999433372
8 21 106250977507865123996
9 44 16793767022807396063419059642469257036652437
10 86 86283684091087378792197424215901018542592662739248420412325877158692469321558575676264
11 179 307954157519416310480198744646044941023074672212201592336666825190665221680585174032224052483643672228833833882969097257874618885560816937172481027939807172551469349507844122611544

精度可以估算為大於最大系數 的二進位制形式的大小


週期
1 1 1 16
2 1 1 16
3 3 2 16
4 6 2 16
5 15 7 16
6 27 13 16
7 63 34 32
8 120 67 64
9 252 144 128
10 495 286 300
11 1023 596 600

以下是 Maxima 程式碼,用於計算

period_Max:11;
/* ----------------- definitions -------------------------------*/
/* 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=0,c] */
GiveG[p]:=
block(
[f:divisors(p),t:1],
g,
f:delete(p,f),
if p=1 
then return(Fn(p,0,c)),
for i in f do t:t*GiveG[i],
g: Fn(p,0,c)/t,  
return(ratsimp(g))
 )$
/* degree of function */
GiveDegree(_function,_var):=hipow(expand(_function),_var);  
log10(x) := log(x)/log(10);
/* ------------------------------*/
file_output_append:true; /* to append new string to file */
grind:true;  
for period:1 thru period_Max step 1 do
(
g[period]:GiveG[period], /* function g */
d[period]:GiveDegree(g[period],c), /* degree of function g */
cf[period]:makelist(coeff(g[period],c,d[period]-i),i,0,d[period]), /* list of coefficients */
cf_max[period]:apply(max, cf[period]), /* max coefficient */
disp(cf_max[period]," ",ceiling(log10(cf_max[period]))),
s:concat("period:",string(period),"  cf_max:",cf_max[period]),
stringout("max_coeff.txt",s)/* save output to file as Maxima expressions */
);


另請參見

尋找中心的圖形方法

[編輯 | 編輯原始碼]

所有這些方法都顯示了週期為 n 及其因子的中心。

動畫顯示顏色與迭代 1-20 的絕對值成正比
第一次迭代的絕對值
顏色顯示第五次迭代落在哪個象限
顏色與 zn 的大小成正比
[編輯 | 編輯原始碼]
  • 顏色與 zn 的大小成正比[18]
  • 引數平面掃描點 c,使得臨界點的軌道消失 [19]
  • YouTube 影片:曼德布洛特振盪 [20]
顏色顯示 zn 落在哪個象限
[編輯 | 編輯原始碼]
九重分解。影像和原始碼

這是曼德布洛特集合外部的徑向第 n 重分解(與 LSM/M 的第 n 重分解進行比較)

使用 4 種顏色,因為有 4 個象限

  • re(z_n) > 0 且 im(z_n) > 0(第一象限)
  • re(z_n) < 0 且 im(z_n) > 0(第二象限)
  • re(z_n) < 0 且 im(z_n) < 0(第三象限)
  • re(z_n) > 0 且 im(z_n) < 0(第四象限)。

“... 當四種顏色在某個地方相遇時,你就有了 q_n(c) 的零點,即週期為 n 的因子的中心。此外,淺色或深色表示 c 是否在 M 中。”(Wolf Jung)

以下是用 C 程式碼計算點 c = cx + cy * i 的 8 位顏色的片段

/* initial value of orbit = critical point Z= 0 */
                       Zx=0.0;
                       Zy=0.0;
                       Zx2=Zx*Zx;
                       Zy2=Zy*Zy;
                       /* the same number of iterations for all points !!! */
                       for (Iteration=0;Iteration<IterationMax; Iteration++)
                       {
                           Zy=2*Zx*Zy + Cy;
                           Zx=Zx2-Zy2 +Cx;
                           Zx2=Zx*Zx;
                           Zy2=Zy*Zy;
                       };
                       /* compute  pixel color (8 bit = 1 byte) */
                       if ((Zx>0) && (Zy>0)) color[0]=0;   /* 1-st quadrant */
                       if ((Zx<0) && (Zy>0)) color[0]=10; /* 2-nd quadrant */
                       if ((Zx<0) && (Zy<0)) color[0]=20; /* 3-rd quadrant */
                       if ((Zx>0) && (Zy<0)) color[0]=30; /* 4-th quadrant */

還可以從 Wolf Jung 的程式 Mandel 原始碼的 mndlbrot.cpp 檔案中的 mndlbrot 類中的 quadrantqn 函式中找到 cpp 程式碼 [21]

與標準迴圈的不同之處是

  • 沒有溢位測試(例如,對於 abs(zn) 的最大值)。這意味著每個點都有相同的迭代次數。對於較大的 n,可能會發生溢位。此外,不能使用測試 Iteration==IterationMax
  • 迭代限制相對較小(例如,從 IterationMax = 8 開始)

另請參閱 Sage 中的程式碼 [22]

使用高斯波函式法 [23]
[編輯 | 編輯原始碼]

牛頓法

[編輯 | 編輯原始碼]

牛頓法尋找一個根

中心列表

[編輯 | 編輯原始碼]
  • James Artherton 發現的週期為 1-20 的中心(總和 = 1 045 239)。使用自定義軟體,採用 320 位定點數學(大約 92 位十進位制數)和網格的四叉樹細分來提高根密度較高的區域的搜尋解析度。這可以適應天線尖端的非常精細的解析度。 [24]
  • Robert P Munafo 發現的最大島嶼 [25] [26]
  • 透過追蹤外部射線找到的所有周期不超過 16 的島嶼的資料庫(週期、島嶼性、帶角度的內部地址、下外部角度分子、分母、上分子、分母、方向、大小、中心實部、虛部) [27]
  • Jay R. Hill 計算的元件的真實中心 [28]
  • 中心列表的集合 [29]
  • 島嶼
  • madore.org/~david math/mandpoints.dat


列表

  • c = 0.339410819995598 -0.050668285162643 i 週期 = 11,扭曲的小曼德布洛特
  • c = 0.325589509550660 -0.038047880934756 i 週期 = 12,扭曲的小曼德布洛特
  • c = 0.314559489984000 -0.029273969079000 i 週期 = 13,扭曲的小曼德布洛特
  • c = 0.272149607027528 +0.005401159465460 i 週期 = 22,扭曲的小曼德布洛特 [30]
85 0.355534 -0.337292
134 -1.7492046334590113301 -2.8684660234660531403e-04 
268 -1.7492046334594190961 -2.8684660260955536656e-04 
272 3.060959246472445584e-01 2.374276727158354376e-02 
204 3.06095924633099439e-01 2.3742767284688944e-02 
204 3.06095924629285095e-01 2.37427672645622342e-02 
136 3.06095924643046857e-01 2.374276727237906e-02 
68 3.06095924629536442e-01 2.37427672749394494e-02 
134 -1.7492046334590114 -2.8684660234659111e-04 
267 -1.1822493706369402373694114253571e-01 6.497492977188836930425943612155e-01
268 -1.182402951560276787014475129283e-01 6.4949165134945441813936036487738e-01
18 -0.814158841137593 0.189802029306573 
32 0.2925755 -0.0149977
52 -0.22817920780250860271129306628202459167994 1.11515676722969926888221122588497247465766

mandelbrot-solver

[編輯 | 編輯原始碼]

這是一個適用於 Ubuntu 的軟體包。一個基於 MPSolve 的曼德布洛特多項式求解器。一個多項式求根器,可以確定具有保證的包含半徑的任意精度近似值。它支援利用多項式結構,如稀疏性,並允許多項式隱式定義或以某些非標準基底定義。該二進位制檔案(作為 MPSolve 軟體包中自定義多項式實現的示例提供)利用了此多項式族體的特殊結構來開發一個高效的求解器,該求解器表現出實驗性的 O(n^2) 複雜度。

此軟體包包含 mpsolve 對曼德布洛特多項式的專門化。

示例

mandelbrot-solver
Usage: mandelbrot-solver n [starting_file]

Parameters: 
 - n is the level of the Mandelbrot polynomial to solve

 - starting_file is an optional file with the approximations that shall be 
                 use as starting points.
a@zelman:~$ mandelbrot-solver 2
 (-1.22561166876654e-01,  7.44861766619744e-01)
 (-1.75487766624669e+00,  0.00000000000000e+00)
 (-1.22561166876654e-01, -7.44861766619744e-01)
a@zelman:~$ mandelbrot-solver 1
 (-1.00000000000000e+00,  0.00000000000000e+00)
a@zelman:~$ mandelbrot-solver 3
 ( 2.82271390766914e-01,  5.30060617578525e-01)
 (-1.56520166833755e-01,  1.03224710892283e+00)
 (-1.94079980652948e+00,  1.26679600613393e-16)
 (-1.00000000000000e+00,  0.00000000000000e+00)
 (-1.31070264133683e+00, -7.22823506473338e-18)
 (-1.56520166833755e-01, -1.03224710892283e+00)
 ( 2.82271390766914e-01, -5.30060617578525e-01)

似乎

 level = period -1

然後結果就是中心。

  1. ^ 人們推測這些多項式(用於中心)是不可約的,但這尚未得到證明 - Wolf Jung

參考文獻

[編輯 | 編輯原始碼]
  1. 核 - 來自 Robert Munafo 的曼德布洛特集合詞彙表和百科全書,(c) 1987-2015。
  2. 複雜動力學中的手術,Carsten Lunde Petersen 線上論文
  3. 牛頓法在實踐中:有效地找到一百萬次多項式的所有根,由 DIERK SCHLEICHER 和 ROBIN STOLL
  4. 維基百科中的除數
  5. FF 主題:曼德布洛特多項式根挑戰(閱讀 5233 次)
  6. Piers Lawrence 和 Rob Corless:曼德布洛特多項式和矩陣。 1; 048; 575 週期為 21 的根
  7. IAP 海報 2013
  8. Robin Stoll 和 Dierk Schleicher 的“牛頓法在實踐中 II:迭代細化牛頓法和尋找某些非常大次數的多項式的所有根的近似最優複雜度”
  9. 維基百科:MPSolve
  10. G Alvarez、M Romera、G Pastor 和 F Montoya 的“透過確定曼德爾布羅特集合的雙曲分量中心”。《混沌、孤子與分形》1998 年,第 9 卷第 12 期,第 1997-2005 頁
  11. Myrberg,P J,“二次實多項式的迭代 III”,芬蘭科學院學報,A. I 編號 348,1964。
  12. 可以計算 100 萬個根的 evaldraw 指令碼,由 knighty 建立
  13. LEONARDO ROBOL 的“多項式和長期方程的求根演算法”
  14. mathoverflow 問題:曼德爾布羅特集合的 Gleason 週期 n 多項式的判別式
  15. 維基百科上的不可約因子
  16. V Dolotin 和 A Morozow:關於基本域的形狀,或者為什麼曼德爾布羅特集合由幾乎理想的圓組成?
  17. maxima 手冊:fpprec
  18. J.D. Long 的“R 中的曼德爾布羅特集合”
  19. Fractal Stream 幫助檔案
  20. youtube 影片:曼德爾布羅特振盪
  21. Wolf Jung 的跨平臺 cpp 程式 Mandel
  22. Christopher Olah 的“逃逸時間分形的形成”
  23. _3__The_Dark Exploding the Dark Heart of Chaos by Chris King March-Dec 2009
  24. James Artherton 的“分形龍”
  25. Robert P. Munafo 的“最大島嶼”,2010 年 10 月 21 日
  26. mrob:largest-islands.txt
  27. Claude 透過跟蹤外部射線找到的所有島嶼的資料庫,這些島嶼的週期最多為 16
  28. Jay R. Hill 計算的分量的實數中心
  29. mset 中心
  30. math.stackexchange 問題:迷你曼德爾布羅特是完全的副本嗎?
華夏公益教科書