GPU支持分形痴迷

我是分形迷。现在我有了一个支持我沉迷的GPU。图形处理器最初是用来加速图形的，但是MATLAB^®使用它们来加速计算。分形¹是需要大量计算的图形。GPU的完美人选。

GPU重新激发了我对Mandelbrot集合本身的兴趣，并促进了三种分形变体的工作:“燃烧的船”，“权力之塔”和牛顿方法的全局收敛行为。

我的GPU是NVIDIA的^®Titan V安装在一个比笔记本电脑大的独立外围外壳中，提供独立的电源和冷却。它计算这些分形的速度大约是CPU的300倍，完全改变了我与我的曼德布洛特程序。我介绍了gpuArray进入程序，现在可以使用与我的屏幕分辨率相当的网格分辨率和迭代，直到图像中的精细细节变得可见。

Mandelbrot集合内部

曼德尔布罗特发生的事就留在曼德尔布罗特。Mandelbrot包含一个简单的复数迭代，从初始点\(z_0\)开始。曼德尔布罗特集是复平面上的一个区域，由值\(z_0\)组成，其中轨迹由

\[z_{k+1} = z_{k}^2 + z_{0}， k= 1,2，…\]

仍然有界。就是这样。这就是整个定义。令人惊讶的是，如此简单的定义可以产生如此迷人的复杂性。

上面显示的彩色图像是曼德尔布罗特集合的鸟瞰图。它没有分辨率来显示在集合边界外的丰富的细节结构的边缘。事实上，Mandelbrot集合有微小的细丝延伸到边缘区域，尽管在这个视图中很难看到边缘。

图1是Mandelbrot集合的几何示意图。最大的部件是一个心形的心形，以参数方程曲线为界

\[z = e^{2it}/ 2 - e^{2it}/4， -\pi \le t \le \pi \]

心线左侧是半径为1 / 4的圆形圆盘。心脏线和椎间盘在红点处接触。稍后会详细介绍那个红点。

围绕着这两个部件的是无限多的近乎圆形的圆盘，圆盘上的圆盘，直径越来越小。最小圆盘的确切位置和形状只能通过计算来确定。

最近已经证明了Mandelbrot集在数学上是连通的，但连通区域有时太薄，以至于我们无法在图形屏幕上解析它，甚至无法在合理的时间内计算它。

计算曼德布洛特

我们与曼德尔布罗特联系在一起的迷人的模式来自于集合之外的边缘区域。

这是计算的核心函数。输入是一个复标量起点z0和一个实标量参数深度．输出是一个计数kz．如果迭代因为z在半径为2的圆盘外，那么z注定成为无限而伯爵kz可作为颜色图的索引。另一方面，如果z幸存了下来深度迭代，然后声明它在Mandelbrot集中。

函数kz = mandelbrot(z0,depth) z = z0;Kz = 0;而(z*conj(z) <= 4) & (kz <= depth) kz = kz + 1;Z = Z * Z + z0;结束结束

让我们在复平面上生成一个点网格，覆盖一个宽度为平方的平面w，以复点为中心佐．

网格= 512;S = w*(-1:1/格:1);[u,v] = meshgrid(s+real(zc)，s+imag(zc));

现在，在CPU上运行的程序和在GPU上运行的程序之间的唯一区别就来了。要在CPU上生成网格，

Z0 = u + i*v;

对于GPU，

z0 = gpuArray(u + i*v);

然后，对于任何一个处理器，语句

Kz = arrayfun(@mandelbrot,z0);

应用标量函数曼德布洛特赋给数组中的所有元素z0．在CPU上，这是一个向量化的double为在网格上循环。在GPU上，语句被分解成数百个并行运行的独立任务。

现在让我们来看两张来自曼德尔布罗特边缘的非凡图像。

边缘地带

这个被称为海马谷的区域位于中心心脏和左侧的椎间盘之间。实际上有两个谷，一个在实轴上方，一个在实轴下方。它们在图1中的红点处相遇。只要稍微发挥一下想象力，您就可以想象图2中所示的小型海洋生物。稍后我们将看到，山谷也包含隐藏的π。

因为曼德尔布罗特集是自相似的，它包含了无限个微型曼德尔布罗特，每个微型曼德尔布罗特的形状都与大曼德尔布罗特的形状相同。图3所示的放大率为10¹⁰．

埋π

1991年，当时还是科罗拉多州立大学研究生的戴夫·波尔发现了一个了不起的结果。波尔当时正在海马谷调查曼德尔布罗特迭代的行为。山谷在接近轴时变得更窄，它们在(-3/4,0)处相遇，即图1中的红点。

我们可以在一个小网格上重复Boll的计算，网格的中心在离轴的-3/4 +处易的极小价值y虚数单位我．

Y = 1.0e-07 zc = -3/4 + Y *i

使网格刚好大到可以接触到轴。

宽度= 2*y网格= 4

选择深度与…成反比y，这让它变得巨大。

深度= 4/y

使用这些参数，运行上面的代码。它产生

Kz = 40000000 40000000 40000000 40000000 40000000 40000000 40000000 40000000 40000000 40000000 40000000 31415926 40000000 40000000 40000000 40000000 40000000 40000000 40000000 20943951 40000000 40000000 40000000 40000000 40000000 15707963 40000000 40000000 40000000

看看网格中心的迭代计数。看到一个熟悉的值了吗?

这不是侥幸。Aaron Klebanoff在2001年的一篇论文《Mandelbrot集合中的π》中分析了类似的计算方法，即在心线前端的尖点处进行计算。

接下来，是Mandelbrot迭代的一个奇怪变体。

燃烧的船

燃烧的船来自一个奇怪的迭代:

\[z_{k+1} = F(z_k) + z_{0}， k= 1,2，…\]

在哪里

\[F(z) = (|Re(z)| + i |Im(z)|)^2\]

我说这很奇怪，因为函数F(z)不是解析函数。我对这个迭代感兴趣，因为下面的图片有不可思议的相似之处。

如图4所示，初始域是一个宽度为3.5的正方形，居中位置为-0.5-0.5我．我插入了一个箭头指向船尾流的目标区域。

将船的尾流放大500倍，达到-1.861- 0.002我．应用骨色彩图使它看起来是冷的，而不是燃烧。图5显示了1915年南极探险家欧内斯特·沙克尔顿(Ernest Shackleton)的船耐力号(Endurance)被冻结在威德尔海(Weddell Sea)冰中的照片旁边的分形结果。

权力之塔

从任何复数开始，\(z\)，反复求它的幂。

\[z, z^z, z^{z^z}， . .\]

我们可以把它表示成一个迭代。从\(y_0 = 1\)开始，让

\[y_{k+1} = z^{y_k}\]

如果\(z\)太大，这个迭代将膨胀到无穷大。对于某个\(z\)，它会收敛到一个有限极限。例如，如果\(z = \sqrt 2\)，则\(y_k\)将收敛于2。最有趣的情况是\(y_k\)接近一个循环。例如，如果\(z\)接近2.5+4\(i\)，则周期长度为7。

2.4684 + 4.0754i -0.6216 + 0.3634i 0.2603 - 0.0184i 1.4868 + 0.3613i -3.4877 + 6.1054i 7.7632e-06 - 2.6617e-06i 1.0000 + 0.0000i 2.4684 + 4.0755i

这个周期长度是“权力之塔”分形的基础。图6显示了整体的分形。

放大10倍⁵并且改变颜色映射会产生如图7所示的图像。

牛顿法

当牛顿方法的起始点不接近函数的零点时，全局行为可能显得不可预测。从复杂平面的起始点区域开始，迭代计算到收敛的等高线图生成发人深省的分形图像。

这种迭代是熟悉的。选择一个函数\(f(z)\)，导数\(f '(z)\)。从\(z_0\)开始，让

\[z_{k+1} = z_{k} - f(z_{k})/f '(z_{k})\]

这最终会收敛到0 (f\)。数一数接近目标所需要的迭代次数。

有许多功能(和颜色地图)可供选择。我最喜欢的三次多项式是

\[f(z) = z^3 - 2z - 5\]

图8显示了将复平面划分为三个区域，其中迭代收敛于三次立方体的三个零之一。在这些区域之间是分形作用强烈的区域，如图中黑色部分所示。

图9显示了牛顿方法求零点的全局行为

\[f(z) = tan(sin(z)) - sin(tan(z))\]

这个函数有无穷多个零和无界的一阶导数。

图10的函数是

\[f(z) = z\;sin(1/z) \]

最突出的蓝色区域围绕着零点±1/π。

这里描述的许多图像对我来说都是新的——我以前从未见过它们。与GPU的交互实验使它们成为可能。

¹分形一词是由波兰/法国/美国数学家Benoit Mandelbrot(1924-2010)创造的，他的大部分职业生涯都在纽约州约克敦高地的IBM沃森研究中心度过