松鼠的窝
首发于松鼠的窝
10561 怎样在球面上「均匀」排列许多点?(下)

10561 怎样在球面上「均匀」排列许多点?(下)

.

  在这篇文章中,我们首先来分析一下公式 (4,5) 产生的平面点阵是怎样受转角\phi影响的,然后把结论拓展到公式 (1~3) 产生的球面点阵上去。

  ★ \begin{align} z_n &= (2n-1)/N - 1 & (1) \\ x_n &= \sqrt{1 - z_n^2} \cdot \cos(2 \pi n \phi) & (2) \\ y_n &= \sqrt{1 - z_n^2} \cdot \sin(2 \pi n \phi) & (3) \end{align}\begin{align} x_n &= \sqrt{n} \cos(2\pi n \phi) & (4)\\ y_n &= \sqrt{n} \sin(2\pi n \phi) & (5) \end{align}

  分析的过程主要依赖连分式(Continued fraction),在此先列举一些连分式的基础知识。任何一个实数a可以展开成如下形式的连分式:

其中各层的分子均为 1;a_0a的整数部分(取下整);其它的各个a_i都是正整数。如果a是有理数,连分式就会像上面一样是有限的;如果a是无理数,连分式则是无限的。上面的式子写起来很麻烦,常常简记为[a_0; a_1, a_2, \ldots, a_n],即写成数列的形式,a_i也称为连分式的项。

  举几个例子:有理数 10/23 写成连分式是0 + \frac{1}{\displaystyle 2 + \frac{1}{\displaystyle 3 + \frac{1}{3}}},是有限的,简记为[0; 2, 3, 3]。而黄金分割比\phi = (\sqrt{5} - 1)/2 \approx 0.618写成连分式是0 + \frac{1}{\displaystyle 1 + \frac{1}{\displaystyle 1 + \frac{1}{\displaystyle 1 + \frac{1}{\ddots}}}},是无限的,简记为[0; 1, 1, 1, \ldots]。连分式的项不一定有规律,例如圆周率\pi的连分式为[3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, \ldots]A001203 - OEIS)。

  把a的连分式在某一项处截断,即可得到a的一个有理逼近,称为丢番图逼近Diophantine approximation)。例如,把\pi的连分式截断为[3;7],则可以得到「约率」22/7;截断为[3; 7, 15, 1],则可以得到「密率」355/113(请读者自行验证)。一般地,若在连分式中保留到a_k这一项,化简后可以得到最简分数p_k/q_k,其中分子和分母可以通过以下的递推式得到:

  ★ \begin{align} p_k &= a_k p_{k-1} + p_{k-2} & (6) \\ q_k &= a_k q_{k-1} + q_{k-2} & (7) \end{align}

递推式的初值为p_0 = a_0, q_0 = 1; p_{-1} = 1; q_{-1} = 0

  丢番图逼近在某种意义下是a最佳逼近。这里「最佳」的意义是:定义p_k/q_k的逼近误差\epsilon_k = |q_ka-p_k|,则不存在另一个分母不超过q_k的分数,逼近误差比\epsilon_k小。如果把a的各个丢番图逼近排成一个数列\frac{p_0}{q_0}, \frac{p_1}{q_1}, \frac{p_2}{q_2}, \ldots,则各项的逼近误差\epsilon_0, \epsilon_1, \epsilon_2, \ldots是递减的。

  (注意,在此我们把误差定义成了\epsilon_k = |q_ka-p_k|,而不是\left| a - \frac{p_k}{q_k} \right|。前一种定义更适合本文的讨论,因为若把a看成每产生一个点旋转的圈数,那么误差\epsilon_k的含义就是:连续产生q_k个点转过的总圈数,与p_k圈差了多少。)

  有了上面的基础知识,我们来分析一些平面点阵。比如,我们取\phi = 0.617。作为准备,我们把 0.617 表示成连分式:[0; 1, 1, 1, 1, 1, 1, 3, 21],并计算出它的各阶丢番图逼近及误差:

  下面是\phi = 0.617时,最靠近原点的 400 个点。可以看到,在原点附近,点阵形成螺旋,共有 13 条旋臂;而在外围,点阵呈放射状,共有 47 条射线。不难发现,13 和 47 都是丢番图逼近的分母。我们把螺旋、射线这种有规律的排列称为「模式」(pattern),并用丢番图逼近来命名,例如本图中的两种模式分别称为 8/13 模式和 29/47 模式。

  200 号点大约处在两种模式的分界线上。观察 200 号点附近的情况:离它最近的点有 153 号、187 号、213 号、247 号,它们的编号与 200 的差恰好也是 13 和 47。把这些点用线连起来,可以发现编号公差为 13 的点列(174 - 187 - 200 - 213 - 226)落在一条 8/13 模式(螺旋)上,而编号公差为 47 的点列(106 -153 - 200 - 247 - 294)落在一条 29/47 模式(射线)上。之所以 8/13 模式产生了明显的旋转,而 29/47 模式不产生明显的旋转,是因为 29/47 的逼近误差远小于 8/13 模式,即:连续产生 47 个点转过的角度与 29 圈的误差,远小于连续产生 13 个点转过的角度与 8 圈的误差。另外还可以发现,在 200 号点附近,8/13 模式上的点的间距随点的编号在增大,而 29/47 模式上的点间距在减小。因此,随着点编号的增大,8/13 模式越来越模糊,而 29/47 模式越来越清晰,导致螺旋被射线所取代。

  有人会说:「140、166、234、260 这几个点,离 200 号点的距离也不远呀!」这几个点的编号与 200 的差是 34 和 60 —— 发现了没有,它们是 13 与 47 的和与差!把它们对应的 21/34 模式和 37/60 模式也画出来可以发现,这些模式上点的间距比 8/13 模式和 29/47 模式大,因此不如后者显著,可以忽略。这正是因为 21/34 和 37/60 不是 0.617 的最佳逼近 —— 21/34 的逼近误差大于 8/13,37/60 的逼近误差大于 29/47。

  经过上面的分析,我们可以得到如下的经验:

    • 每个最佳逼近p_k/q_k,对应着一种比较显著的模式;
    • 模式上的点间距越小,模式越显著;
    • 同一个模式上的点间距随位置会发生变化,导致最显著模式的切换。

为了发掘模式切换的规律,下面我们来定量地计算模式上点的间距。

  考虑n号点,以及它所处的p_k/q_k模式。在该模式上,下一个点的编号是n+q_k。把这两个点的编号代入公式 (4,5),就能得到两点的坐标和距离了。但是这样得到的公式太复杂,我们换一个角度来思考。n号点距原点的距离为r=\sqrt{n},在n号点附近,编号每增加 1,距原点的距离增加\frac{\mathrm{d}r}{\mathrm{d}n} = \frac{1}{2\sqrt{n}} = \frac{1}{2r}。当q_k \ll n时,从n号点到n+q_k号点,径向距离就是\frac{q_k}{2r}。而从n号点到n+q_k号点转过的总圈数为q_k\phi,它大约是p_k整圈,误差为\epsilon_k = |q_k \phi - p_k|,于是这两点间的切向距离就是2\pi r \epsilon_k。当\epsilon_k \ll 1时,可以把这个切向距离看成与径向垂直的直线段,由此得到p_k/q_k模式在距原点r处的点间距:

  ★ \begin{align} d_k &= \sqrt{\left(\frac{q_k}{2r}\right)^2 + (2\pi r \epsilon_k)^2} & (8) \end{align}

这里我选择距原点距离r而不是点的编号n为自变量,显得更直观。

  如果选定某一种模式(即固定q_k, \epsilon_k),则d_kr的函数。根号下的第一项随r递减,第二项随r递增,整体上则呈现先减后增的「V」字形。d_k的最小值为d_{k,\min} = \sqrt{2 \cdot \frac{q_k}{2r} \cdot 2\pi r \epsilon_k} = \sqrt{2\pi q_k \epsilon_k},它在\frac{q_k}{2r} = 2\pi r \epsilon_kr_{k,\min} = \sqrt{\frac{q_k}{4\pi\epsilon_k}}处取得。由此可知,从原点向外,各种模式会先变得越来越显著,这是因为半径r的增长越来越慢,使得模式上的点靠得越来越近;然后会变得越来越模糊,因为转角q_k\phip_k整圈的微小偏差被越来越大的半径放大成越来越远的距离。这个规律也可以通俗地描述成:每个模式先是形成越来越清晰的射线,然后射线弯曲成螺旋,最后「被风吹散」。根据r_{k,\min} = \sqrt{\frac{q_k}{4\pi\epsilon_k}}k值越大的模式,分母q_k越大,而误差\epsilon_k越小,故使得该模式最显著的半径r_{k,\min}越大。也就是说,k值越大的模式,就会在离原点越远处显现,这造成了模式的交替。

  我们用上面\phi = 0.617时的平面点阵检验一下上面的结论。下表计算了各个模式最显著处的半径r_{k,\min}(称为「最显著半径」)和此处的点间距d_{k,\min}(称为「最密点间距」)。k \le 5的各个模式(包括 0/1, 1/1, 1/2, 2/3, 3/5, 5/8),误差都太大,同时r_{k,\min}太小、d_{k,\min}太大,所以在图上显现不出来。8/13 模式在半径约为 7 处(50 号点附近)最显著,这与图示相符;29/47 模式在半径约为 61 处最显著,但前文的图中只画出了半径为 20 以内的部分,所以这一点就看不出来了。从表中可以看到 8/13 模式的最密点间距为 1.310,而 29/47 模式的最密点间距为 0.543,所以若把点阵继续画下去, 29/47 模式会变得比 8/13 模式更显著。最后还有一个模式 617/1000,它将在 29/47 模式消散后显现出来。并且因为 617/1000 就是\phi的精确值,此模式将永不消散,并越来越显著。

  令d_k = d_{k+1},可以求出模式p_k/q_kp_{k+1}/q_{k+1}过渡处的半径。仍然研究\phi = 0.617的例子,取k=6,可以求得模式 8/13 与 29/47 之间的过渡发生在半径为 13.090 处(170 号点附近)。我们前面观察过的 200 号点离这儿也不远,从图上可以看到这里确实发生了模式的过渡。同样,取k=7,可以求得模式 29/47 到 617/1000 的过渡发生在半径为 281.939 处,也就是说,要画到 80000 号点,才能观察到这个过渡。

  为了让大家一饱眼福,我把点阵画到了 30 万个点,并截取了横轴正半轴周围的局部。请点开大图观察,否则图片上会产生摩尔纹(Moiré pattern)。从图中可以清楚地看到原点附近的 8/13 模式(螺旋),它在半径为 13 处过渡为 29/47 模式,后者在半径 60 左右最显著,且比 8/13 模式更显著。29/47 模式的竖纹在半径 250 处基本消散,到半径 350 处开始产生 617/1000 模式的横纹。

  当\phi = 0.617时各模式的点间距随半径的变化可以总结成下图,注意横、纵轴均为对数刻度。从图中可以清楚地看出 8/13、29/47、617/1000 三种主要模式的最显著半径、最密点间距,以及模式之间过渡的位置。一图胜千言,下文中我将主要使用这种图(称为模式图)来说明模式的变化。

  下面,我们来探究连分式中项的大小对模式的影响。我们选取\phi = \pi - 3 \approx 0.1416。这个值的连分式为[0; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, \ldots],其中 15、292 这样的大项,会带给我们新的发现。

  先画个点阵看看:

图中可以观察到两种模式:中心的模式有 7 条腿,外围的模式腿太多了,数不清。

  按上文给出的公式,计算\phi的丢番图逼近序列、各个逼近的误差,以及各个模式的最显著半径和最密点间距,并作出模式图,如下所示。

  从模式图中观察到,1/7 和 16/113 这两个模式十分显著,它们在半径为 30 处左右发生过渡,这与点阵图完全相符。注意,在这两个模式发生过渡的时候,有一个 15/106 模式被完全掩盖过去了。16/113 这个模式会一直持续到半径接近 10000 的时候才会过渡到 4703/33215 这个模式,哪怕画大图也无缘查看了。

  1/7 和 16/113 这两个模式,正是连分式中 15 和 292 这两个大项产生的。这里的「大」,其实并不是指绝对数值很大,而是指相对于周围的项来说很大。回忆模式最密点间距公式d_{k,\min} = \sqrt{2\pi q_k \epsilon_k},其中q_k是递增的,\epsilon_k是递减的。再回忆丢番图逼近中分母的递推式q_k = a_k q_{k-1} + q_{k-2}:当a_k是一个大项时,q_kq_{k-1}相比,就会猛然增大。但当a_k很大时,若把连分式截断至a_{k-1},那么舍弃的\frac{1}{a_k + \ddots}就很小,所以\epsilon_{k-1}与之前的\epsilon_{k-2}相比就会猛然减小。误差的减小比分母的增大早一步发生,这就使得d_{k-1, \min}明显小于相邻模式的最密点间距,产生一个十分显著的模式p_{k-1}/q_{k-1}。另外观察模式图可以发现,各个模式的「V」字形的两臂斜率都是一样的,所以越是显著的模式,覆盖的半径范围往往也越大。

  上面的计算告诉我们,连分式中每出现一个大项,在它之前截断就会产生一个十分显著的模式。那么,若要不产生显著的模式,连分式中就不能出现大项,所以令所有的a_k都等于 1 就是最好的选择了。这就给出了最优的转角 —— 黄金角。当\phi = (\sqrt{5} - 1) / 2 \approx 0.618时,模式图如下(注意,各模式的分子、分母都是菲波那契数):

  可以看出,各个模式的显著程度都差不多,并且模式切换频繁。观察下面点阵图的中心,似乎能看出许多层层叠叠的花瓣,这正是模式频繁切换的体现。到了外围,模式持续的长度就长了些(注意模式图的横坐标是对数刻度),能看出一些模式了;这些模式都呈螺旋而不是射线状,因为各阶丢番图逼近的误差都不小。

  模式切换频繁正是点阵混乱和密集的原因:混乱是因为任何一种模式(规律)都不会持续很久;密集是因为在任一半径处,都并没有单独一种模式的点间距比别的模式小得多,所以不会出现大的缝隙。

  可以证明,对于黄金分割比\phi,丢番图逼近p_k/q_k的误差\epsilon_k \approx \frac{1}{\sqrt{5}q_k},于是各模式的最密点间距就都约等于\sqrt{2 \pi / \sqrt{5}} \approx 1.676。这其实是用形如 (4,5) 的公式能生成的点阵中,模式最密点间距的「上界」—— 加引号是因为,可能出现个别模式的最密点间距大于这个值,但这样的个别模式至多只有有限个。

  话说回来,公式 (4,5) 并不是平面上点阵密铺问题的最优解。在保证同样点阵密度(每\pi面积上有一个点)的情况下,正方形网格的最密点间距是d_\text{square} = \sqrt{\pi} \approx 1.772,而最优解六边形网格的最密点间距是d_\text{hex} = \sqrt{2\pi/\sqrt{3}} \approx 1.905(如下图)。这两种网格满足均匀、密集,但不满足混乱,因为同一种排列规律遍布了整个平面。从 1.905 到 1.676,就是「混乱」需要付出的代价。

  现在来看一下我在试探过程中发现的另一个不错的\phi值:\phi = \sqrt{2} - 1 \approx 0.414。此时的点阵图与模式图如下:

  把\phi = \sqrt{2} - 1 \approx 0.414写成连分式,是[0; 2, 2, 2, \ldots],各项都相等,所以没有「大项」。模式图上也显示出,各个模式的最密点间距都相等,且模式切换频繁。但是,与黄金分割比相比,\phi = \sqrt{2} - 1 \approx 0.414的各个模式的最密点间距偏小(约为\sqrt{\pi/\sqrt{2}} \approx 1.490),故模式还是稍微明显了一些;模式之间的切换也不如黄金分割比频繁(注意两张模式图的横、纵坐标范围是相同的)。黄金分割比以微弱的优势战胜\sqrt{2} - 1

  上面我们详细地分析了\phi = 0.617,\,\, \pi - 3,\,\, \frac{\sqrt{5} - 1}{2},\,\, \sqrt{2} - 1各值时的点阵图与模式图,得到了如下一些结论:

    • 每个点转过的圈数\phi的每个丢番图逼近对应着一种模式;
    • 随半径增加,每种模式先变得显著,然后消散,存在一个最显著半径
    • 分母越大的模式,最显著半径也越大,所以各模式会依次显现
    • 当连分式中有大项的时候,大项之前的那个丢番图逼近对应的模式会非常显著,且持续的长度也很长;
    • 黄金分割比的连分式中所有项都是 1,这导致没有显著的模式模式切换频繁,于是点阵显得密集而混乱。

  终于到了把结论拓展到球面的时候了。在平面上,我们选择了到原点的距离r作为自变量;在球面上,我们选择纬度\theta作为自变量。设整个点阵的点数为N,考虑模式p_k/q_k上的n号点和n+q_k号点。在纬度为\theta处,这两个点在经线方向上的距离为\frac{2q_k}{N\cos\theta},在纬线方向上的距离为2 \pi \epsilon_k \cos\theta。由此可得p_k/q_k模式在竖坐标为z处的点间距:

  ★ \begin{align} d_k &= \sqrt{\left(\frac{2q_k}{N \cos\theta}\right)^2 + (2\pi \epsilon_k \cos\theta)^2} & (9) \end{align}

此式的最小值为d_{k,\min} = \sqrt{8\pi q_k \epsilon_k / N},在\theta_{k,\min} = \pm \arccos\sqrt{\frac{q_k}{N\pi\epsilon_k}}处取得。作为对比,平面上p_k/q_k模式的最小点间距为d_{k,\min} = \sqrt{2\pi q_k \epsilon_k},在r_{k,\min} = \sqrt{\frac{q_k}{4\pi\epsilon_k}}处取得。两个最小点间距的公式是相似的,只是系数不同;但最显著位置的公式则有形式上的差别:球面上的公式多了个反余弦。反余弦函数是单调递减的,且定义域有限,这说明,各个模式是从两极到赤道依次显现的,而有些高阶的模式,哪怕到了赤道处也来不及显现。除此以外,平面点阵的各个结论都适用于球面点阵。

  例如,下面是\phi = 0.617时,在球面上取 1000 个点生成的点阵,以及相应的模式图。可以看到 8/13 模式(螺旋)在纬度 64 度时最显著,并在 35 度左右过渡成 29/47 模式(射线)。但在低纬度区域,8/13 模式并没有消散得很厉害,它与 29/47 模式共同组成了方形网格。617/1000 模式来不及显现出来。

  下图是\phi取黄金分割比时,由 1000 个点组成的球面点阵和模式图。与平面的情况类似,这里没有特别显著的模式,且模式切换频繁,所以点阵的分布密集且混乱。不过,由于 arccos 定义域的影响,21/34 这个模式从北纬 40 度一直到南纬 40 度都是最显著的;但与此同时,相邻的 13/21 模式和 34/55 模式的显著程度也差不多,所以在低纬度区域的点阵中能找到正方形或六边形网格的影子。

编辑于 2017-09-27

文章被以下专栏收录