分形/复平面迭代/动态外射线

图像

• 落在alpha不动点上的射线
• 
  双可达
• 
  三可达

• 落在临界点和临界值上的射线
• 
  1/4
• 
  1/6
• 
  9/56
• 
  129/16256

理论

令$\Psi _{c}\,$ 为从闭合单位圆盘${\overline {\mathbb {D} }}$ 的补集（外部）到填充的Julia集$\ K_{c}$ 的补集的映射（保角映射，同构）。

$\Psi _{c}:{\hat {\mathbb {C} }}\setminus {\overline {\mathbb {D} }}\to {\hat {\mathbb {C} }}\setminus K_{c}$

其中${\hat {\mathbb {C} }}$ 表示黎曼球面（扩展复平面）。

令$\Psi _{c}=\Phi _{c}^{-1}\,$ 表示

• 逆Boettcher映射
• 黎曼映射

角度为$\theta \,$ 的外射线记为${\mathcal {R}}_{\theta }^{K}$ 是

• 在$\Psi _{c}\,$ 下的直线${\mathcal {R}}_{\theta }=\{\left(r\cdot e^{2\pi i\theta }\right):\ r>1\}$ 的像

${\mathcal {R}}_{\theta }^{K}=\Psi _{c}({\mathcal {R}}_{\theta })$

• 具有相同外部角 $\theta$ 的填充后的 Julia 集外部的点集

${\mathcal {R}}_{\theta }^{K}=\{z\in {\hat {\mathbb {C} }}\setminus K_{c}:\arg(\Phi _{c}(z))=\theta \}$

重要问题

哪些是独特的射线？

• "落在 Julia 集的关键点（或包含关键点的球体的根点）上的射线 t 对该 Julia 集是唯一的，因此是识别 Julia 集的一种方式。因此，我用 t 来命名 Julia 集 Jc。在余下部分中，我只能处理 t=p/q 为有理数的 Julia 集。"
• 临界前周期多项式通常由 落在临界值的外部射线 的角度 $\theta$ 进行参数化

如何逼近外部射线？

• 用于 Julia 集的二进制分解方法 (BDM/J)
• 场线 的标量场（势）
• SAM = 条带平均方法 基本上是计算角度的一种廉价方式。

• 
• 
  BDM 的边界
• 
  SAM

在 BDM 图像中，角度（以圈数测量）的外部射线

$angle=(k/2^{n})~~{\mbox{mod }}~1\,$

可以看作子集的边界。

当两个外部动态射线落在同一点时

问题：

• Q1：给定 Julia 集中的一个点，有多少条射线落在该点上？
• Q2：给定两条射线 $R_{\theta _{1}}$ 和 $R_{\theta _{2}}$，在什么条件下它们落在同一点？
  • 当外部射线落在同一点时
  • 着陆等效性 = 哪些射线落在同一点？
  • 数值标准 = 它们具有相同的（数值）着陆点
• Q3：哪些射线支持相同的 Fatou 分量？

定理 1.2 设 f 为度数 d ≥ 2 的多项式，其 Julia 集局部连通。如果两个角度 $\theta _{1}$ 和 $\theta _{2}$ 关于临界肖像具有相同的行程，那么外部射线 $R_{\theta _{1}}$ 和 $R_{\theta _{2}}$ 的着陆点要么重合，要么属于 Fatou 域的边界，该边界最终被迭代到 Siegel 盘上。

如何找到落在临界值 z= c 上的外部射线的角度呢？

这取决于具体情况。

• 例如，当参数 c 是一个 Misiurewicz 点时，可以使用组合方法。
• 大多数情况下，你首先构造角度，然后确定参数。
• 但是，如果你想从 Julia 集或 Hubbard 树中读取角度，你可以按照 Mandel 程序 中的 Demo 4 第 1 和 2 页的说明进行。

构造 填充 Julia 集 的脊柱

构造脊柱的算法由 A. Douady 描述

• 连接 $-\beta \,$ 和 $\beta \,$,
• (待完成)

• 

绘制外部射线是什么意思呢？

这意味着

• 计算（近似）射线的某些点
• 用线段连接这些点

这将给出射线 ${\mathcal {R}}$ 的近似值。

绘制 动态外部射线

算法

• 牛顿法
• 逆向迭代
• Boettcher 映射的逆映射
• 追踪角度为 t 的外部射线（曲线）

逆向迭代

对复二次多项式进行反向迭代，并适当选择前像

周期为 2-40 的落在 α 固定点上的周期性外部射线。使用反向迭代生成。

这种方法已经被许多人使用，并被 Thierry Bousch 证明。

Wolf Jung 的 C++ 代码可以在 qmnplane.cpp 文件中的 QmnPlane::backray() 函数中找到（参见 Mandel 程序 的源代码）。

• 周期性角度的射线（最简单的情况）

我们将通过一个例子来解释它

首先选择外部角度 $t\,$ （以圈数为单位）。周期性射线的外部角度是一个有理数。

$t={\frac {1}{3}}\,$

计算外部角度在 倍增映射 下的周期。

因为 "1/3 的倍增结果是 2/3，而 2/3 的倍增结果是 4/3，它与 1/3 同余"

$1\equiv 2{\pmod {3}}.\,$

或者

${\frac {2}{3}}\equiv {\frac {1}{3}}{\pmod {1}}.\,$

因此外部角 $t={\frac {1}{3}}\,$ 在倍增映射下周期为2。

从无穷远处的两个点开始（在共轭平面上）

在射线 1/3 上有一个点 $w=R*e^{2\pi i/3}\,$

在射线 2/3 上有一个点 $w=R*e^{4\pi i/3}\,$.

在无穷远附近 $z=w\,$ 因此可以切换到动力学平面（博特切尔共轭）

向后迭代（从两个可能性中适当选择） 在射线 1/3 上的点会到射线 2/3，再回到 1/3，依此类推。

在 C 语言中是 

/* choose one of 2 roots: zNm1 or -zNm1  where zN = sqrt(zN - c )  */
if (creal(zNm1)*creal(zN) + cimag(zNm1)*cimag(zN) <= 0) zNm1=-zNm1;

或者在 Maxima CAS 

if (z1m1.z01>0) then z11:z1m1 else z11:-z1m1;

必须将点集分成 2 个子集（2 条射线）。绘制这两个集合中的一个，并将这些点连接起来。它将近似于射线。

• 对于预周期角的射线（待完成）

/*

compute last point ~ landing point
of the dynamic ray for periodic angles ( in turns )

 gcc r.c -lm -Wall -march=native

landing point of ray for angle = 1 / 15 = 0.0666666666666667 is = (0.0346251977103306 ;  0.4580500411138030 ) ; iDistnace = 18 
landing point of ray for angle = 2 / 15 = 0.1333333333333333 is = (0.0413880816505388 ;  0.5317194187688231 ) ; iDistnace = 17 
landing point of ray for angle = 4 / 15 = 0.2666666666666667 is = (-0.0310118081927549 ;  0.5440125864026020 ) ; iDistnace = 17 
landing point of ray for angle = 8 / 15 = 0.5333333333333333 is = (-0.0449867688014234 ;  0.4662592852362425 ) ; iDistnace = 18

*/

// https://gitlab.com/adammajewski/ray-backward-iteration
#include <stdio.h>
#include <stdlib.h> // malloc
#include <math.h> // M_PI; needs -lm also 
#include <complex.h>

/* --------------------------------- global variables and consts ------------------------------------------------------------ */
#define iPeriodChild 4 // iPeriodChild of secondary component joined by root point

// - --------------------- functions ------------------------

/* 
   principal square  root of complex number 
   wikipedia  Square_root

   z1= I;
   z2 = root(z1);
   printf("zx  = %f \n", creal(z2));
   printf("zy  = %f \n", cimag(z2));
*/
double complex root(double x, double y)
{ 
  
  double u;
  double v;
  double r = sqrt(x*x + y*y); 
  
  v = sqrt(0.5*(r - x));
  if (y < 0) v = -v; 
  u = sqrt(0.5*(r + x));
  return u + v*I;
}

// This function only works for periodic  angles.
// You must know the iPeriodChild n before calling this function.
// draws all "iPeriodChild" external rays 
// commons  File:Backward_Iteration.svg
// based on the code by Wolf Jung from program Mandel
// http://www.mndynamics.com/

int ComputeRays( //unsigned char A[],
			    int n, //iPeriodChild of ray's angle under doubling map
			    int iterMax,
                            double Cx, 
                            double Cy,
                            double dAlfaX,
                            double dAlfaY,
                            double PixelWidth,
                            complex double zz[iPeriodChild] // output array

			    )
{
  double xNew; // new point of the ray
  double yNew;
  
  const double R = 10000; // very big radius = near infinity
  int j; // number of ray 
  int iter; // index of backward iteration

  double t,t0; // external angle in turns 
  double num, den; // t = num / den
  
  double complex zPrev;
  double u,v; // zPrev = u+v*I
 
  int iDistance ; // dDistance/PixelWidth = distance to fixed in pixels

  
  /* dynamic 1D arrays for coordinates ( x, y) of points with the same R on preperiodic and periodic rays  */
  double *RayXs, *RayYs;
  int iLength = n+2; // length of arrays ?? why +2

  //  creates arrays :  RayXs and RayYs  and checks if it was done
  RayXs = malloc( iLength * sizeof(double) );
  RayYs = malloc( iLength * sizeof(double) );
  if (RayXs == NULL || RayYs==NULL)
    {
      fprintf(stderr,"Could not allocate memory");
      getchar(); 
      return 1; // error
    }

   // external angle of the first ray 
   num = 1.0;
   den = pow(2.0,n) -1.0;
   t0 = num/den; // http://fraktal.republika.pl/mset_external_ray_m.html
   t=t0;
   // printf(" angle t = %.0f / %.0f = %f in turns \n", num, den, t0);

  //  starting points on preperiodic and periodic rays 
  //  with angles t, 2t, 4t...  and the same radius R
  for (j = 0; j < n; j++)
    { // z= R*exp(2*Pi*t)
      RayXs[j] = R*cos((2*M_PI)*t); 
      RayYs[j] = R*sin((2*M_PI)*t);
      //
      // printf(" %d angle t = = %.0f / %.0f =  %.16f in turns \n", j, num , den,  t); 
      //
      num *= 2.0;
      t *= 2.0; // t = 2*t
      if (t > 1.0) t--; // t = t modulo 1
      
    }
  //zNext = RayXs[0] + RayYs[0] *I;

  // printf("RayXs[0]  = %f \n", RayXs[0]);
  // printf("RayYs[0]  = %f \n", RayYs[0]);

  // z[k] is n-periodic. So it can be defined here explicitly as well.
  RayXs[n] = RayXs[0]; 
  RayYs[n] = RayYs[0];

  //   backward iteration of each point z
  for (iter = -10; iter <= iterMax; iter++)
    { 
     	
      for (j = 0; j < n; j++) // n +preperiod
	{ // u+v*i = sqrt(z-c)   backward iteration in fc plane 
	  zPrev = root(RayXs[j+1] - Cx , RayYs[j+1] - Cy ); // , u, v
	  u=creal(zPrev);
	  v=cimag(zPrev);
                
	  // choose one of 2 roots: u+v*i or -u-v*i
	  if (u*RayXs[j] + v*RayYs[j] > 0) 
	    { xNew = u; yNew = v; } // u+v*i
	  else { xNew = -u; yNew = -v; } // -u-v*i

	  // draw part of the ray = line from zPrev to zNew
	 // dDrawLine(A, RayXs[j], RayYs[j], xNew, yNew, j, 255);
                
	  
	  //  
	  RayXs[j] = xNew; RayYs[j] = yNew;
                
	

                
	} // for j ...

          //RayYs[n+k] cannot be constructed as a preimage of RayYs[n+k+1]
      RayXs[n] = RayXs[0]; 
      RayYs[n] = RayYs[0];
          
      // convert to pixel coordinates 
      //  if z  is in window then draw a line from (I,K) to (u,v) = part of ray 
   
      // printf("for iter = %d cabs(z) = %f \n", iter, cabs(RayXs[0] + RayYs[0]*I));
     
    }

  // Approximate end of ray by straight line to it's landing point here = alfa fixed point
 // for (j = 0; j < n + 1; j++)
  //  dDrawLine(A, RayXs[j],RayYs[j], dAlfaX, dAlfaY,j, 255 );

  // this check can be done only from inside this function
  t=t0;
  num = 1.0;
  for (j = 0; j < n ; j++)
    {

      zz[j] = RayXs[j] +  RayYs[j] * I; // save to the output array
      // Approximate end of ray by straight line to it's landing point here = alfa fixed point
      //dDrawLine(RayXs[j],RayYs[j], creal(alfa), cimag(alfa), 0, data); 
      iDistance = (int) round(sqrt((RayXs[j]-dAlfaX)*(RayXs[j]-dAlfaX) +  (RayYs[j]-dAlfaY)*(RayYs[j]-dAlfaY))/PixelWidth);
      printf("last point of the ray for angle = %.0f / %.0f = %.16f is = (%.16f ;  %.16f ) ; Distance to fixed = %d  pixels \n",num, den, t, RayXs[j], RayYs[j],  iDistance);
      num *= 2.0;
      t *= 2.0; // t = 2*t
      if (t > 1) t--; // t = t modulo 1
    } // end of the check

  // free memmory
  free(RayXs);
  free(RayYs);

  return  0; //  
}

int main()
{
  
  complex double l[iPeriodChild];
  int i;
  // external angle in turns = num/den; 
  double num = 1.0;
  double den = pow(2.0, iPeriodChild) -1.0;

   ComputeRays( iPeriodChild, 
                   10000, 
                   0.25, 0.5, 
                   0.00, 0.5,
                   0.003,
                   l ) ;

  
  printf("\n see what is in the output array : \n");

  for (i = 0; i < iPeriodChild ; i++) {
       printf("last point of the ray for angle = %.0f / %.0f = %.16f is = (%.16f ;  %.16f ) \n",num, den, num/den, creal(l[i]), cimag(l[i]));
       num *= 2.0;}

  return 0;
}

运行它 

./a.out

输出 

last point of the ray for angle = 1 / 15 = 0.0666666666666667 is = (0.0346251977103306 ;  0.4580500411138030 ) ; Distance to fixed = 18  pixels 
last point of the ray for angle = 2 / 15 = 0.1333333333333333 is = (0.0413880816505388 ;  0.5317194187688231 ) ; Distance to fixed = 17  pixels 
last point of the ray for angle = 4 / 15 = 0.2666666666666667 is = (-0.0310118081927549 ;  0.5440125864026020 ) ; Distance to fixed = 18  pixels 
last point of the ray for angle = 8 / 15 = 0.5333333333333333 is = (-0.0449867688014234 ;  0.4662592852362425 ) ; Distance to fixed = 19  pixels

 see what is in the output array : 
last point of the ray for angle = 1 / 15 = 0.0666666666666667 is = (0.0346251977103306 ;  0.4580500411138030 ) 
last point of the ray for angle = 2 / 15 = 0.1333333333333333 is = (0.0413880816505388 ;  0.5317194187688231 ) 
last point of the ray for angle = 4 / 15 = 0.2666666666666667 is = (-0.0310118081927549 ;  0.5440125864026020 ) 
last point of the ray for angle = 8 / 15 = 0.5333333333333333 is = (-0.0449867688014234 ;  0.4662592852362425 )

射线上的点正在向后移动

• 当它离茱莉亚集合很远时，速度很快
• 靠近茱莉亚集合时速度很慢（经过 50 次迭代后，点之间的距离 = 0 像素）

参见计算示例（这里 像素大小 = 0.003） 

# iteration distance_between_points_in_pixels
0 	 3300001 
1 	 30007 
2 	 2296 
3 	 487 
4 	 179 
5 	 92 
6 	 54 
7 	 34 
8 	 23 
9 	 18 
10 	 14 
11 	 11 
12 	 9 
13 	 7 
14 	 6 
15 	 5 
16 	 5 
17 	 4 
18 	 4 
19 	 3 
20 	 3 
21 	 3 
22 	 3 
23 	 2 
24 	 2 
25 	 2 
26 	 2 
27 	 2 
28 	 2 
29 	 2 
30 	 1 
31 	 1 
32 	 1 
33 	 1 
34 	 1 
35 	 1 
36 	 1 
37 	 1 
38 	 1 
39 	 1 
40 	 1 
41 	 1 
42 	 1 
43 	 1 
44 	 1 
45 	 1 
46 	 1 
47 	 1 
48 	 1 
49 	 1 
50 	 1 
51 	 1 
52 	 1 
53 	 1 
54 	 1 
55 	 1 
56 	 0 
57 	 0 
58 	 0 
59 	 0 
60 	 0 

可以选择像素大小不同的点 

#iteration    distance(z1,z2)   distance (z,alfa)
       0 	 3300001 	 33368
       1 	   30007 	  3364
       2 	    2296 	  1074
       3 	     487 	   591
       4 	     179 	   413
       5 	      92 	   321
       6 	      54 	   267
       7 	      34 	   234
       8 	      23 	   211
       9 	      18 	   193
      10 	      14 	   179
      11 	      11 	   169
      12 	       9 	   160
      13 	       7 	   153
      14 	       6 	   146
      15 	       5 	   141
      16 	       5 	   136
      17 	       4 	   132
      18 	       4 	   128
      19 	       3 	   125
      20 	       3 	   122
      21 	       3 	   119
      22 	       3 	   117
      23 	       2 	   115
      24 	       2 	   112
      25 	       2 	   110
      26 	       2 	   109
      27 	       2 	   107
      28 	       2 	   105
      29 	       2 	   104
      30 	       1 	   102
      31 	       1 	   101
      32 	       1 	   100
      33 	       1 	    99
      34 	       1 	    97
      35 	       1 	    96
      36 	       1 	    95
      38 	       2 	    93
      40 	       2 	    92
      42 	       2 	    90
      44 	       2 	    88
      46 	       1 	    87
      48 	       1 	    86
      50 	       1 	    84
      52 	       1 	    83
      54 	       1 	    82
      56 	       1 	    81
      59 	       1 	    80
      62 	       1 	    78
      65 	       1 	    77
      68 	       1 	    76
      71 	       1 	    75
      74 	       1 	    74
      78 	       1 	    73
      82 	       1 	    71
      86 	       1 	    70
      90 	       1 	    69
      95 	       1 	    68
     100 	       1 	    67
     105 	       1 	    66
     111 	       1 	    65
     117 	       1 	    64
     124 	       1 	    63
     131 	       1 	    62
     139 	       1 	    60
     147 	       1 	    59
     156 	       1 	    58
     166 	       1 	    57
     177 	       1 	    56
     189 	       1 	    55
     202 	       1 	    54
     216 	       1 	    53
     231 	       1 	    52
     247 	       1 	    51
     265 	       1 	    50
     285 	       1 	    49
     307 	       1 	    48
     331 	       1 	    47
     358 	       1 	    46
     388 	       1 	    45
     421 	       1 	    44
     458 	       1 	    43
     499 	       1 	    42
     545 	       1 	    41
     597 	       1 	    40
     655 	       1 	    39
     721 	       1 	    38
     796 	       1 	    37
     881 	       1 	    36
     978 	       1 	    35
    1090 	       1 	    34
    1219 	       1 	    33
    1368 	       1 	    32
    1542 	       1 	    31
    1746 	       1 	    30
    1986 	       1 	    29
    2270 	       1 	    28
    2608 	       1 	    27
    3013 	       1 	    26
    3502 	       1 	    25
    4098 	       1 	    24
    4830 	       1 	    23
    5737 	       1 	    22
    6873 	       1 	    21
    8312 	       1 	    20
   10157 	       1 	    19
   12555 	       1 	    18
   15719 	       1 	    17
   19967 	       1 	    16
   25780 	       1 	    15
   33911 	       1 	    14
   45574 	       1 	    13
   62798 	       1 	    12
   89119 	       1 	    11
  131011 	       1 	    10
  201051 	       1 	     9
  325498 	       1 	     8
  564342 	       1 	     7
 1071481 	       1 	     6
 2308074 	       1 	     5
 5996970 	       1 	     4
21202243 	       1 	     3
136998728 	       1 	     2

可以观察到，从像素 12 移动到 alfa 附近的像素 11 需要 27 000 次迭代。计算到 alfa 附近 1 像素内的点需要：2m1.236 秒

使用柯蒂斯·麦克穆伦的逆博特切尔映射绘制动态外部射线

使用麦克穆伦方法绘制带有外部射线的茱莉亚集合

该方法基于柯蒂斯·麦克穆伦的 C 程序 及马特亚兹·艾拉特的 Pascal 版本

这里有两个平面 

• w-平面（或 f0 平面）
• z-平面（fc 平面的动力学平面）

该方法包含 3 个主要步骤 

• 计算圆形外部射线在角 $\theta \,$ 和不同半径（光栅化）上的某些 w 点

$w_{0}=R*e^{2\pi it}\,$ 其中 $1<R<\infty \,$

• 使用逆博特切尔映射将 w 点映射到 z 点 $\Psi _{c}=\Phi _{c}^{-1}\,$

$z_{0}=\Psi _{c}(w_{0})=f_{c}^{-n}(w_{0}^{2^{n}})\,$

• 绘制 z 点（并使用线段连接它们（线段是直线的一部分，由两个不同的端点界定））

第一步和最后一步很容易，但第二步不容易，需要更多解释。

光栅化

对于给定的外部射线，在 $f_{0}\,$ 平面 中，射线的每个点都有 

• 常数 $t\,$ （外角以转为单位）
• 变量半径 $R\,$

所以 $w\,$ 射线的点由半径 $R\,$ 参数化，可以使用 复数的指数形式 计算。

$w=F_{t}(R)=R*e^{2*\Pi *i*t}\,$

可以使用 **线性尺度** 沿着射线前进。

t:1/3; /* example value */
R_Max:4;
R_Min:1.1;
for R:R_Max step -0.5 thru R_Min do w:R*exp(2*%pi*%i*t);
/* Maxima allows non-integer values in for statement */

这会得到一些 w 点，这些点之间的距离相等。

另一种方法是使用 **非线性尺度**。

为了做到这一点，我们引入浮点 指数 $r\,$，使得

$R=2^{r}\,$

和

$w=G_{t}(r)=2^{r}*e^{2*\Pi *i*t}\,$

为了计算外射线在 $f_{0}\,$ 平面上的角度为 $t\,$ 的一些 w 点，使用以下 Maxima 代码

t:1/3; /* external angle in turns */
/* range for computing  R ; as r tends to 0 R tends to 1 */
rMax:2; /* so Rmax=2^2=4 /
rMin:0.1; /* rMin > 0   */
caution:0.93; /* positive number < 1 ; r:r*caution  gives smaller r */
r:rMax;
unless r<rMin do
( 
 r:r*caution, /* new smaller r */ 
 R:2^r, /* new smaller R */
 w:R*exp(2*%pi*%i*t) /* new point w in f0 plane */
);

在这个方法中，点之间的距离并不相等，而是与填充的 Julia 集边界距离 成反比。

这是好的，因为在这里射线具有更大的 曲率，因此曲线会更加平滑。

映射

从 $f_{0}\,$-平面 到 $f_{c}\,$-平面 的映射包括 4 个小步骤

• 在 $f_{0}\,$ 平面中进行向前迭代

$w_{n}=w_{0}^{2^{n}}=R^{2^{n}}*e^{2\pi i2^{n}t}\,$

直到 $w_{n}\,$ 接近无穷大

$abs(w_{n})>R_{limit}\,$

• 切换平面（从 $f_{0}\,$ 到 $f_{c}\,$)

$z_{n}=w_{n}\,$

（因为这里，接近无穷大：$w=z\,$）

• 在 $f_{c}\,$ 平面上进行反向迭代相同（$n\,$）次迭代次数
• 最后一个点 $z_{0}\,$ 在我们的外部射线上

步骤 1、2 和 4 很容易。第三个不容易。

$z_{n-1}=f_{c}^{-1}(z_{n})={\sqrt {z_{n}-c}}\,$

反向迭代使用复数的平方根。它是一个双值函数，因此反向迭代会产生二叉树。

在没有额外信息的情况下，无法在这样的树中选择正确的路径。为了解决这个问题，我们将使用两件事

• 无穷大吸引域的等度连续性
• $f_{0}\,$ 和 $f_{c}\,$ 平面之间的共轭关系

无穷大吸引域的等度连续性

无穷大吸引域（填充的 Julia 集的补集）包含在正向迭代下趋于无穷大的所有点。

$A_{f_{c}}(\infty )\ {\overset {\underset {\mathrm {def} }{}}{=}}\ \{z\in \mathbb {C} :f_{c}^{(k)}(z)\to \infty \ as\ k\to \infty \}.$

无穷大是一个超吸引不动点，所有点的轨道都具有类似的行为。换句话说，如果两个点的初始位置很接近，则它们的轨道被认为会保持接近。

这就是等度连续性（与正规性比较）。

在 $f_{c}\,$ 平面上，可以使用射线先前点的正向轨道来计算下一个点的反向轨道。

算法的详细版本

• 计算射线的第一个点（从无穷大附近开始，朝 Julia 集移动）

$w=2^{r}*e^{2\Pi it}\,$ 其中 $r=r_{max}\,$

这里可以轻松切换平面

$z=w\,$

这是射线的第一個 z 点。

• 计算射线的下一个 z 点
  • 计算射线的下一个 w 点，对于 $r=r*caution\,$
  • 计算两个点的正向迭代：前一个 z 点和当前 w 点。保存 z 轨迹和最后一个 w 点
  • 切换平面并使用最后一个 w 点作为起点：$z_{n}=w_{last}\,$
  • 使用前一个 z 点的正向轨迹，对新的 $z_{n}\,$ 进行反向迭代，直到新的 $z_{0}\,$
  • $z_{0}\,$ 是射线的下一个 z 点
• 依此类推（下一个点），直到 $r<r_{Min}\,$

Maxima CAS 源代码

/* gives a list of z-points of external ray for angle t in turns and coefficient c */
GiveRay(t,c):=
block(
 [r],
 /* range for drawing  R=2^r ; as r tends to 0 R tends to 1 */
 rMin:1E-20, /* 1E-4;  rMin > 0  ; if rMin=0 then program has infinity loop !!!!! */
 rMax:2,  
 caution:0.9330329915368074, /* r:r*caution ; it gives smaller r */
 /* upper limit for iteration */
 R_max:300,
 /* */
 zz:[], /* array for z points of ray in fc plane */
 /*  some w-points of external ray in f0 plane  */
 r:rMax,
 while 2^r<R_max do r:2*r, /* find point w on ray near infinity (R>=R_max) in f0 plane */
 R:2^r,
 w:rectform(ev(R*exp(2*%pi*%i*t))),
 z:w, /* near infinity z=w */
 zz:cons(z,zz), 
 unless r<rMin do
 ( /* new smaller R */
  r:r*caution,  
  R:2^r,
  /* */
  w:rectform(ev(R*exp(2*%pi*%i*t))),
  /* */
  last_z:z,
  z:Psi_n(r,t,last_z,R_max), /* z=Psi_n(w) */
  zz:cons(z,zz)
 ),
 return(zz)
)$

另请参见

• 参数外部射线
• 外部射线精度 - 重要性

参考文献

• 外部射线的豪斯多夫极限：拓扑图像，作者：Carsten Lunde Petersen 和 Saeed Zakeri