分形/複平面迭代/離散拉格朗日描述符

連續時間動力系統（拉格朗日描述流的方法） 是分析相空間結構的一種方法。這裡 該方法擴充套件到離散動力系統：複平面上的開放對映。

影像

• 週期 1
• 
  花椰菜朱利亞集 - 週期 1 - 拋物線情況
• 
  拋物線情況：場線

• 週期 2
• 
  週期 2 吸引情況
• 
  巴西利卡朱利亞集 - 週期 2 超吸引情況
• 
  胖巴西利卡：週期 2 拋物線 - 場線

• 不連線
• 
  內爆花椰菜

完整原始碼在公共頁面上（點選影像）

關鍵詞

複數

$z=x+yi$

復對映

對映

$z_{n+1}=f(z_{n})=f^{(n+1)}(z_{0})\quad \forall \ n\ \in \mathbb {N} \cup \left\{0\right\}$

黎曼球面

黎曼球面上的點 $s$

$s=(\xi _{1},\xi _{2},\xi _{3})$

立體投影逆變換

從北極到球體下方平面的立體投影的 3D 圖示

立體投影逆變換將複平面的點 $z$ 對映到黎曼球面上的點 $s$

$S^{-1}(z)=\left({\frac {2x}{1+x^{2}+y^{2}}},{\frac {2y}{1+x^{2}+y^{2}}},{\frac {x^{2}+y^{2}-1}{1+x^{2}+y^{2}}}\right)$

所以

$s=S^{-1}(z)=(\xi _{1},\xi _{2},\xi _{3})$

和

 $\xi _{1}={\frac {2x}{1+x^{2}+y^{2}}}$
 
$\xi _{2}={\frac {2y}{1+x^{2}+y^{2}}}$

$\xi _{3}={\frac {x^{2}+y^{2}-1}{1+x^{2}+y^{2}}}$

p-範數

向量 $\mathbf {x} =(x_{1},\ldots ,x_{n})$ 的 $p$-範數（也稱為 $\ell _{p}$-範數）是

$\left\|\mathbf {x} \right\|_{p}:=\sum _{i=1}^{n}\left|x_{i}\right|^{p}$

其中

• 數字 p 是實數 $p\in (0,1]$。它被稱為冪。它會影響奇點（如朱利亞集等分形特徵）附近的梯度的陡峭程度

p-範數用於度量對映 f 在黎曼球面上連續迭代之間的距離

離散拉格朗日描述符 = DLD

  The simple idea is to compute the p-norm version of Lagrangian descriptors, not for the points on the complex plane, but for their projections on the Riemann sphere in the extended complex plane.
  ... in the complex mappings that we consider in this work, the functions that define the dynamics are not invertible, and therefore we will only keep the forward part of the definition.  

DLD

• 是一個標量值
• 累積軌道的 p-範數（= 包含有關軌道歷史的資訊），因此揭示了朱利亞集內部和外部的結構

求和

 summing is what the original paper does ( pauldelbrot)

$D_{p}(z_{0},N)=\sum _{k=0}^{N-1}\left\|s_{k+1}-s_{k}\right\|_{p}=\sum _{j=1}^{3}\sum _{k=0}^{N-1}\left|\xi _{j}^{(k+1)}-\xi _{j}^{(k)}\right|^{p}$

其中

• N 是固定迭代次數
• $z_{0}=x_{0}+y_{0}i\$ 是複平面上有界子集 D 上選取的任何初始條件
• $s$ 是黎曼球面上的一個點：$s=(\xi _{1},\xi _{2},\xi _{3})$

平均

 "averaging keeps the coloring stable if maxiters is changed but can lead to low variation over the image" pauldelbrot

${\frac {D_{p}(z_{0},N)}{N}}$

步驟

對於複平面上的每個點 z

• 計算 DLD（標量值）
• 顏色與 DLD 成比例

子步驟：計算 z 的 DLD

• 在對映 f 下迭代點 z = 計算 zn
• 將每個點 zn 從複平面對映到黎曼球面（反立體投影）
• 對於每個 zn 計算被加數
• ... (待辦事項)

程式碼

UltraFractal

DLD {
; Based on https://arxiv.org/pdf/2001.08937.pdf
; ucl file for UltraFractal by pauldelbrot
init:
  float sum = 0.0
  float lastx = 0.0
  float lasty = 0.0
  float lastz = 0.0
  int i = 0
loop:
  float d = |#z|
  float dd = 1/(d + 1)
  ; Riemann sphere coordinates = (xx, yy,zz)
  float xx = 2*real(#z)*dd
  float yy = 2*imag(#z)*dd
  float zz = (d - 1)*dd     
  :
  IF (i > 0)
    sum = sum + (xx - lastx)^@power + (yy - lasty)^@power + (zz - lastz)^@power
  ENDIF
  i = i + 1
  lastx = xx
  lasty = yy
  lastz = zz
final:
  #index = sum/(i - 1)
default:
  title = "Discrete Langrangian Descriptors"
  param power
    caption = "Power"
    default = 0.25
    hint = "Affects the steepness of the gradient near singularities (fractal features like a Julia set)"
    min = 0.0
  endparam
}

  " Here's the latest version I've been using in UF. It handles escaping points with no-bail formulae via the isInf/isNaN test (puts the Riemann sphere point at the north pole for those), allows averaging or summing (summing is what the original paper does, whereas averaging keeps the coloring stable if maxiters is changed but can lead to low variation over the image), and can use or not use absolute values on the differences being summed." pauldelbrot

DLD {
; Based on https://arxiv.org/pdf/2001.08937.pdf
; ucl file for UltraFractal by pauldelbrot
init:
  float sum = 0.0
  float lastx = 0.0
  float lasty = 0.0
  float lastz = 0.0
  int i = 0
loop:
  float d = |#z|
  float dd = 1/(d + 1)
  float xx = 2*real(#z)*dd
  float yy = 2*imag(#z)*dd
  float zz = (d - 1)*dd     ; Riemann sphere coordinates
  IF (isInf(d) || isNaN(d))
    ; Infinity, or thereabouts
    xx = 0
    yy = 0
    zz = 1
  ENDIF
  IF (i > 0)
    IF(@qabs)
      sum = sum + abs(xx - lastx)^@power + abs(yy - lasty)^@power + abs(zz - lastz)^@power
    ELSE
      sum = sum + (xx - lastx)^@power + (yy - lasty)^@power + (zz - lastz)^@power
    ENDIF
  ENDIF
  i = i + 1
  lastx = xx
  lasty = yy
  lastz = zz
final:
  IF(@qsum)
    #index = sum
  ELSE
    #index = sum/(i - 1)
  ENDIF
default:
  title = "Discrete Langrangian Descriptors"
  param power
    caption = "Power"
    default = 0.25
    hint = "Affects the steepness of the gradient near singularities (fractal features like a Julia set)"
    min = 0.0
  endparam
  param qsum
    caption = "Sum"
    default = false
    hint = "Averages if false, sums if true."
  endparam
  param qabs
    caption = "Abs differences"
    default = false
  endparam
}

Fragmentarium

基於 UF 程式碼的程式碼，由 3Dickulus 為 GLSL 修改和最佳化

#include "Complex.frag"
#include "MathUtils.frag"
#include "Progressive2D.frag"
#info Unveiling Fractal Structure with Lagrangian Descriptors
#info https://fractalforums.org/fractal-mathematics-and-new-theories/28/unveiling-the-fractal-structure-of-julia-sets-with-lagrangian-descriptors/3376/msg20446#msg20446

#group Lagrangian

// Number of iterations
uniform int  Iterations; slider[1,200,1000]
uniform vec3 RGB; slider[(0,0,0),(0.0,0.4,0.7),(1,1,1)]
uniform bool Julia; checkbox[false]
uniform vec2 JuliaXY; slider[(-2,-2),(-0.6,1.3),(2,2)]
uniform float p; slider[0,.6,1]

/* partial pnorm
   input: z, c, p
   output ppn
*/

float ppnorm( vec2 z, vec2 c, float p){

	vec3 s0,s1; // for 2 points on the Riemann sphere
	float d; // denominator
	float ds;

	// map from complex plane to riemann sphere
	// z
	d = z.x*z.x + z.y*z.y + 1.0;
	s0 = vec3(2.0*z,(d-2.0))/d;
	// zn
	d = c.x*c.x + c.y*c.y + 1.0;
	s1 = vec3(2.0*c,(d-2.0))/d;
	// sum
	vec3 ss = pow(abs(s1 - s0),vec3(p));
   ds = ss.x+ss.y+ss.z;

	return ds;
}

// DLD = Discret Lagrangian Descriptior
float lagrangian( vec2 z, vec2 c, float p ){

	int i; // number of iteration
	float d = 0.0; // DLD = sum

	for (i=0; i<Iterations; ++i){
		d += ppnorm(z, c, p); // sum z
		z = cMul(z,z) +c; // complex iteration
		if (cAbs(z) > 1e19 ) break; // exterior : upper limit of float type
	}

	d /= float(i); // averaging not summation

	return d;
}

vec3 color(vec2 c) {
	vec2 z = Julia ? c : vec2(0.,0.);
	if(Julia) c = JuliaXY;
	float co = lagrangian( z, c, p );
	return .5+.5*cos(6.2831*co+RGB);
}

#preset Default
Center = -0.724636541,0.025224931
Zoom = 0.64613535
EnableTransform = false
RotateAngle = 0
StretchAngle = 0
StretchAmount = 0
Gamma = 2.2
ToneMapping = 1
Exposure = 1
Brightness = 1
Contrast = 1
Saturation = 1
AARange = 2
AAExp = 1
GaussianAA = true
Iterations = 20
RGB = 0,0.4,0.7
p = 0.1444322
Julia = false
JuliaXY = -1.05204872,0
Bailout = 1000
#endpreset

#preset Basilica
Center = -0.025346913,-0.013859176
Zoom = 0.561856826
EnableTransform = false
RotateAngle = 0
StretchAngle = 0
StretchAmount = 0
Gamma = 2.2
ToneMapping = 1
Exposure = 1
Brightness = 1
Contrast = 1
Saturation = 1
AARange = 2
AAExp = 1
GaussianAA = true
Iterations = 20
RGB = 0,0.4,0.7
p = 0.1444322
Julia = true
JuliaXY = -1.05204872,0
Bailout = 1000
#endpreset

c

/* partial pnorm 
   input: z , zn = f(z), p
   output ppn

*/
double ppnorm( complex double z, complex double zn, double p){

	double s[2][3]; // array for 2 points on the Riemann sphere
	int j; 
	double d; // denominator 
	double x; 
	double y;
	
	double ds;
	double ppn = 0.0;
	
	// map from complex plane to riemann sphere
	// z
	x = creal(z);
	y = cimag(z);
	d = x*x + y*y + 1.0;
	
	s[0][0] = (2.0*x)/d;
	s[0][1] = (2.0*y)/d;  
	s[0][2] = (d-2.0)/d; // (x^2 + y^2 - 1)/d
	
	// zn
	x = creal(zn);
	y = cimag(zn);
	d = x*x + y*y + 1.0;
	s[1][0] = (2.0*x)/d;
	s[1][1] = (2.0*y)/d;  
	s[1][2] = (d-2.0)/d; // (x^2 + y^2 - 1)/d
	
	// sum 
	for (j=0; j <3; ++j){
		ds = fabs(s[1][j] - s[0][j]);
		ppn += pow(ds,p); // |ds|^p
		}
	return ppn;
}

// DLD = Discret Lagrangian Descriptior
double lagrangian( complex double z0, complex double c, int iMax, double p ){

	int i; // number of iteration
	double d = 0.0; // DLD = sum
	double ppn; // partial pnorm
	complex double z = z0;
	complex double zn; // next z
	
	
	if (cabs(z) < AR || cabs(z +1)< AR) return 5.0; // for z= 0.0 d = inf
	
	
	for (i=0; i<iMax; ++i){
	
        zn = z*z +c; // complex iteration
		ppn = ppnorm(z, zn, p);
		d += ppn; // sum
		//
		z = zn; 
		
		if (cabs(z) > ER ) break; // exterior : big values produces NAN error in ppnorm computing 
		if (cabs(z) < AR || cabs(z +1)< AR) 
			{ // interior
				d = -d;
				break; 
				
			}
		}
	 
	d =  d/((double)i); // averaging not summation
	if (d<0.0) {// interior
		d = 2.5 - d;
	}
	return d; 
}

unsigned char ComputeColorOfDLD(complex double z){

 	
  	int iColor;
  	double d;

  	d = lagrangian(z,c, N,p);
  	iColor = (int)(d*255)  % 255; // color is proportional to d
  
  
  return (unsigned char) iColor;
}

論文

• 用拉格朗日描述符揭示朱利亞集的分形結構，作者：Víctor J. García-Garrido
  • 在《非線性科學與數值模擬通訊》第 91 卷，2020 年 12 月，105417 中，只有一篇匿名評論，但需要付費才能檢視
  • 透過 arxiv.org，免費，但未經同行評審的預印本

www

• fractalforums.org : unveiling-the-fractal-structure-of-julia-sets-with-lagrangian-descriptors/3376/msg19803

參考文獻

V. J. García-Garrido。用拉格朗日描述符揭示朱利亞集的分形結構。《非線性科學與數值模擬通訊》第 91 卷 (2020) 105417。