分形/复平面迭代/离散拉格朗日描述符

连续时间动力系统（拉格朗日描述流的方法） 是分析相空间结构的一种方法。这里 该方法扩展到离散动力系统：复平面上的开放映射。

图像

• 周期 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。