跳轉到內容

分形/複平面迭代/曼德勃羅集內部

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

本書展示瞭如何為繪製引數平面[1](曼德勃羅集[2])編寫不同的演算法,用於復二次多項式。[3]

人們可以在引數平面上找到不同型別的點/集。[4]

此頁面介紹了曼德勃羅集的內部點。[5]

曼德勃羅集內部 - 雙曲分量

[編輯 | 編輯原始碼]

“捕獲時間”演算法:收斂所需的迭代次數

[編輯 | 編輯原始碼]
The “capture-time algorithm” is a natural counterpart for points inside the set to the “escape-time algorithm”. Given some desired tolerance, the orbit P is generated for each point c ∈ C until some point in the orbit is closer than to some previous point in the orbit. The number of iterations needed for this to occur is mapped to a color and displayed at the pixel corresponding to c. Adam Cunningham[6]


李雅普諾夫指數

[編輯 | 編輯原始碼]

數學方程 :[7]

其中

表示 f 相對於 z 的一階導數

另請參閱

  • janthor 的影像和描述[8]
  • Anders Sandberg 的影像[9]


JPBotelho 的 HLSL 程式碼[10]

Shader "Fractals/Coloring Techniques/Escape-Time"
{
	Properties
	{
		_MainTex ("Texture", 2D) = "white" {}

		_Iter ("Iterations", Range(0, 250)) = 100
		_Dividend ("Dividend", Range (0, 0.5)) = 15
		
		_Zoom ("Zoom", Range (0.1, 1000)) = 0.65
		_Position ("Offset", Vector) = (0.4, 0, 0, 0)

		_Background ("Background", Color) = (0, 0.25, 1, 0)
		_Origin ("Origin", Color) = (0, 0, 0, 0)
	}	
	SubShader
	{
		Cull Off ZWrite Off ZTest Always

		Pass
		{
			CGPROGRAM

			#pragma vertex vert
			#pragma fragment frag
			
			#include "UnityCG.cginc"
			#include "Complex.cginc"
			#include "FractalOperations.cginc"

			struct appdata
			{
				float4 vertex : POSITION;
				float2 uv : TEXCOORD0;
			};

			struct v2f
			{
				float2 uv : TEXCOORD0;
				float4 vertex : SV_POSITION;
			};

			sampler2D _MainTex;

			int _Iter;
			fixed _Zoom;
			fixed _Dividend;
			float2 _Position;
			fixed4 _Background;
			fixed4 _Origin;

			v2f vert (appdata v)
			{
				v2f o;
				o.vertex = UnityObjectToClipPos(v.vertex);
				o.uv = v.uv;
				return o;
			}

			fixed4 frag (v2f i) : SV_Target
			{
				float x0 = (ClampScaleX(i.uv) + _Position.x) / _Zoom;
				float y0 = (ClampScaleY(i.uv) + _Position.y) / _Zoom;
				
				float2 z = float2(x0, y0);
				float2 c = float2(x0, y0);
				
				int iteration = 0;

				float l = 0;

				while (IsBounded(z, 40) && iteration < _Iter)
				{	
					l += log (cabs(2 * z));

					z = cmul(z, z);
					z += c;					

					iteration++;
				}

				l /= iteration;

				if (l > 0)
					return _Background;

                float3 color = tanh(l >= 0 ? 
											float3(0, 0.7 * log(1 + l), log(1 + l)) : 
											3 * float3(_Origin.x-l, _Origin.y-l * 0.1, _Origin.z));

				return float4(color + _Dividend, 1);


			}		    
			ENDCG
		}
	}
	CustomEditor "FractalEditor"
}

內部距離估計

[編輯 | 編輯原始碼]
內部距離估計

DEM/M - 方法描述

軌道的絕對值

[編輯 | 編輯原始碼]
# Hypercomputing the Mandelbrot Set? by Petrus H. Potgieter February 1, 2008
n=1000; # For an nxn grid
m=50; # Number of iterations
c=meshgrid(linspace(-2,2,n))\ # Set up grid
+i*meshgrid(linspace(2,-2,n));
x=zeros(n,n); # Initial value on grid
for i=1:m
x=x.^2+c; # Iterate the mapping
endfor
imagesc(min(abs(x),2.1)) # Plot monochrome, absolute
# value of 2.1 is escape

內部水平集

[編輯 | 編輯原始碼]

點的顏色

  • 與 z 在最後一次迭代時的值成正比。
  • 顯示了週期性吸引子的內部水平集。

bof60 的影像位於“分形的美麗”這本書的第 60 頁。方法描述在 bof 的第 63 頁。它僅用於曼德勃羅集的內部點。

點的顏色與

  • 其軌道到原點[11][12]的最小距離成正比
  • 迭代過程中 z 獲得的最小值[13]
  • 闡明瞭原點(臨界點)的迭代在集合內部到原點的最近距離
  • “每個影片幀的每個畫素都代表一個特定的複數 c = a + ib。對於每個連續幀 n,z(c,n) := z(c, n-1)^2 + c 的幅度顯示為這些點 c 中的每個點的灰度強度值:幅度較大的點更白,幅度較小的點更暗。當 n 從 1 增加到 256 時,曼德勃羅集外部的點會迅速飽和為純白色,而曼德勃羅集內部的點會在較暗的強度之間振盪。” 布萊恩·高沃特[14]

距離的水平集是具有相同距離的點集。[15]

if (Iteration==IterationMax)
 /* interior of Mandelbrot set = color is proportional to modulus of last iteration */
 else { /* exterior of Mandelbrot set = black */
  color[0]=0;
  color[1]=0;
  color[2]=0;                           
 }
  • 程式碼片段 : 來自 Gnofract4d 的 fractint.cfrm[16]
bof60 {
 init:
       float mag_of_closest_point = 1e100
 loop:
       float zmag = |z|
       if zmag < mag_of_closest_point
               mag_of_closest_point = zmag
       endif
 final:
       #index = sqrt(mag_of_closest_point) * 75.0/256.0
}


另請參閱

bof61 或原子域

[編輯 | 編輯原始碼]

完整描述

雙曲分量的週期

[編輯 | 編輯原始碼]
雙曲分量的週期

曼德勃羅集的雙曲分量的週期是臨界軌道的極限集的週期。

用於計算週期的演算法

  • 直接從動力學平面上的臨界點 z = 0.0 的迭代中檢測週期
  • “快速而骯髒”演算法:檢查是否,則將 c 點用顏色 n 著色。這裡 n 是吸引軌道的週期,eps 是吸引點周圍圓的半徑 = 數值計算的精度
  • “如果正確實現,基於區間算術的方法能夠找到相當大的 n 的所有周期 n 迴圈。” (ZBIGNIEW GALIAS)[17]
  • Floyd 的迴圈查詢演算法[18]
  • 蜘蛛演算法
  • 原子域,BOF61
  • 週期檢測


內部檢測

[edit | edit source]

如果以下所有內容都成立,則畫素很有可能是內部[19]

  • 畫素被標記為內部(黑色)
  • 所有周圍的畫素都被標記為內部(黑色)
  • 所有黑色畫素具有相同的週期

內部座標和乘數對映

[edit | edit source]
使用乘數對映計算的曼德布羅集的元件
曼德布羅集 - 乘數對映

定義

克勞德·海蘭德-艾倫的演算法

  • 檢查 c
    • 當 c 在曼德布羅集之外時
      • 現在放棄
      • 或使用外部座標
    • 當 c 不在外部(在內部或邊界上)時:對於每個週期 p,從 1 開始並遞增
      • 找到週期點 z0,使得 fp(z0,c)=z0,使用牛頓法在一個復變數中
      • 透過在 z0 處對 fp 關於 z 的一階導數進行求值來找到 b
      • 如果 |b|≤1,則返回 b,否則繼續下一個 p

計算

[edit | edit source]

對於週期:[23]

  • 1 到 3 可以使用顯式方程[24]
  • >3 必須使用數值方法來找到

週期 1

[edit | edit source]

邊界方程開始

 c+(w/2)^2-w/2=0;

並針對 w 求解

(%i1) eq1:c+(w/2)^2-w/2=0;
                                                                                                              2
                                                                                                             w    w
(%o1)                                                                                                        -- - - + c = 0
                                                                                                             4    2
(%i2) solve(eq1,w);
(%o2)                                                                                        [w = 1 - sqrt(1 - 4 c), w = sqrt(1 - 4 c) + 1]
(%i3) s:solve(eq1,w);
(%o3)                                                                                        [w = 1 - sqrt(1 - 4 c), w = sqrt(1 - 4 c) + 1]
(%i4) s:map(rhs,s);
(%o4)                                                                                            [1 - sqrt(1 - 4 c), sqrt(1 - 4 c) + 1]

所以

 w = w(c) =  1.0 - csqrt(1.0-4.0*c)

週期 2

[edit | edit source]
 w = 4.0*c + 4;

週期 3

[edit | edit source]
 

它可以使用 Maxima CAS 求解

(%i1) e1:c^3 + 2*c^2 - (w/8-1)*c + (w/8-1)^2 = 0;

                      3      2        w       w     2
(%o1)                c  + 2 c  + (1 - -) c + (- - 1)  = 0
                                      8       8
(%i2) solve(e1,w);
(%o2) [w = (- 4 sqrt((- 4 c) - 7) c) + 4 c + 8, w = 4 sqrt((- 4 c) - 7) c + 4 c + 8]

數值逼近

[edit | edit source]
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;
}


另請參閱

內角

[edit | edit source]
用徑向角著色的曼德布羅集的內部

Renato Fonseca 的方法:[25] “集合中的一個點 c 的色調等於引數

(適當縮放,以便最終得到 0 到 255 範圍內的數字)。數字 z_nmax 是 z 序列中計算的最後一個數字。”

另請參閱


Fractint

[edit | edit source]

Fractint:顏色引數:INSIDE=ATAN

顏色透過確定最後一次迭代值相對於實軸的角度(以度為單位)並使用絕對值來確定。此功能應與 periodicity=0[26] 一起使用

內部射線

[edit | edit source]

沿著內部射線從雙曲引數到拋物線引數[27]


變化而 保持不變時, 沿著 內部射線 移動。 [28] 它被用作曼德勃羅集內部的 路徑


 
double complex Give_c(double t, double r, int p)
{
	/*
	input:
	InternalRadius = r in [0,1] 
  	InternalAngleInTurns = t in range [0,1]
  	p = period
  	
  	output = c = complex point of 2D parameter plane  
  	*/
  	

	complex double w = 0.0;
	complex double c = 0.0;
	
	t = t*2*M_PI; // from turns to radians
  	// point of unit circle
  	w = r* cexp(I*t);
  		
	// map circle to component
	switch (p){
	
	case 1: c = (2.0*w - w*w)/4.0; break;
	case 2: c = (w -4.0)/ 4.0; break;
  
	}
	return c; 
}


/* find c in component of Mandelbrot set 
 uses complex type so #include <complex.h> and -lm 
 uses code by Wolf Jung from program Mandel
 see function mndlbrot::bifurcate from mandelbrot.cpp
 http://www.mndynamics.com/indexp.html

  */
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int period)
{
  //0 <= InternalRay<= 1
  //0 <= InternalAngleInTurns <=1
  double t = InternalAngleInTurns *2*M_PI; // from turns to radians
  double R2 = InternalRadius * InternalRadius;
  double Cx, Cy; /* C = Cx+Cy*i */
  switch ( period ) {
    case 1: // main cardioid
      Cx = (cos(t)*InternalRadius)/2-(cos(2*t)*R2)/4; 
      Cy = (sin(t)*InternalRadius)/2-(sin(2*t)*R2)/4; 
      break;
   case 2: // only one component 
      Cx = InternalRadius * 0.25*cos(t) - 1.0;
      Cy = InternalRadius * 0.25*sin(t); 
      break;
  // for each period  there are 2^(period-1) roots. 
  default: // safe values
      Cx = 0.0;
      Cy = 0.0; 
    break; }

  return Cx+ Cy*I;
}

// draws points to memory array data
int DrawInternalRay(double InternalAngleInTurns, unsigned int period, int iMax, unsigned char data[])
{

   complex double c;
   double InternalRadius;
   double RadiusStep; // between radius of points 
   int i; // number of point to draw
      
  RadiusStep = 1.0/iMax;
   
  for(i=0;i<=iMax;++i){ 
   InternalRadius = i * RadiusStep;
   c = GiveC(InternalAngleInTurns, InternalRadius, period);
   DrawPoint(c,data);
  }

return 0;
}

示例:角度為 1/6 的主心形的內部射線。

內角

射線的半徑

單位圓的內部半徑點

將點 對映到引數平面

對於 ,這是主心形的方程。

內部曲線

[編輯 | 編輯原始碼]

保持不變而 變化時, 沿著內部曲線移動。

/* find c in component of Mandelbrot set 
 uses complex type so #include <complex.h> and -lm 
 uses code by Wolf Jung from program Mandel
 see function mndlbrot::bifurcate from mandelbrot.cpp
 http://www.mndynamics.com/indexp.html
*/
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int period)
{
  //0 <= InternalRay<= 1
  //0 <= InternalAngleInTurns <=1
  double t = InternalAngleInTurns *2*M_PI; // from turns to radians
  double R2 = InternalRadius * InternalRadius;
  double Cx, Cy; /* C = Cx+Cy*i */
  switch ( period ) {
    case 1: // main cardioid
      Cx = (cos(t)*InternalRadius)/2-(cos(2*t)*R2)/4; 
      Cy = (sin(t)*InternalRadius)/2-(sin(2*t)*R2)/4; 
      break;
    case 2: // only one component 
      Cx = InternalRadius * 0.25*cos(t) - 1.0;
      Cy = InternalRadius * 0.25*sin(t); 
      break;
    // for each period  there are 2^(period-1) roots. 
    default: // safe values
      Cx = 0.0;
      Cy = 0.0; 
    break;
  }

  return Cx+ Cy*I;
}

// draws points to memory array data
int DrawInternalCurve(double InternalRadius , unsigned int period,  int iMax, unsigned char data[])
{
  complex double c;
  double InternalAngle; // in turns = from 0.0 to 1.0
  double AngleStep;
  int i;
  // int iMax =100;
   
  AngleStep = 1.0/iMax;
   
  for (i=0; i<=iMax; ++i) { 
    InternalAngle = i * AngleStep;
    c = GiveC(InternalAngle, InternalRadius, period);
    DrawPoint(c,data);
  }

  return 0;
}

分量的中心

[編輯 | 編輯原始碼]

更多教程和程式碼

[編輯 | 編輯原始碼]


教程

參考資料

[編輯 | 編輯原始碼]
  1. 維基百科中的引數平面
  2. 維基百科中的 Mandelbrot 集
  3. 維基百科中的復二次多項式
  4. reenigne 部落格:mandelbrot-set-taxonomy
  5. 透過 A Cunningham 展示 Mandelbrot 集的內部結構(帶有 python 3 程式和程式碼)
  6. Adam Cunningham 展示 Mandelbrot 集的內部結構
  7. 2013 年 10 月 4 日,Didier Gonze 的邏輯方程
  8. janthor 的 Lyapunov 指數和 Mandelbrot 集
  9. Anders Sandberg 的影像
  10. github 庫 JPBotelho:Fractal-Megacollection(Unity 的 HLSL 著色器)
  11. Fractint:各種選項和演算法
  12. Java™ Number Cruncher:Java 程式設計師的數值計算指南,作者 Ronald Mak
  13. Firefly 應用程式幫助,作者 Terry W. Gintz
  14. Mandelbrot 振盪,作者 Brian Gawalt
  15. Fractint 文件,作者 Noel Giffin
  16. gnofract4d
  17. 使用區間方法對電子電路中的週期軌道進行嚴格研究,作者 Zbigniew Galias
  18. Milan 的 Mandelbrot 集繪製
  19. fractalforums.org:使用級數逼近確定跳過的最佳迭代次數
  20. Mandelbrot 集內部座標,作者 Claude Heiland-Allen
  21. 內部距離渲染實踐,作者 Claude Heiland-Allen
  22. math.stackexchange 問題:測試 Mandelbrot 球中週期為 n 的成員資格/1151953#1151953
  23. Brown 方法,作者 Robert P. Munafo,2003 年 9 月 22 日。
  24. 精確座標,作者 Robert P. Munafo,2003 年 9 月 22 日。
  25. Mandelbrot 集,作者 Renato Fonseca
  26. fractint 顏色引數
  27. 沿著內射線的雙曲引數到拋物線引數,作者 Yi-Chiuan Chen 和 Tomoki Kawahira
  28. 維基百科中的內部射線
  29. ASCII 圖形
華夏公益教科書