绘制分形

分形是一种美丽且引人注目的数学形式，通常可以由一套相对简单的指令创建。在自然界，分形可以在不同的位置找到，例如海岸线、贝壳和蕨类植物，甚至被用于创建某些类型的卫星天线。分形的数学概念已被了解了一段时间，但直到 1970 年代才开始受到真正的重视，因为当时计算机图形学取得了进步，并且一些偶然的发现让像 Benoît Mandelbrot 这样的研究人员偶然发现了分形所具有的真正令人惊叹的可视化效果。

今天，我们将学习如何绘制这些美丽的可视化效果，并且将在我们熟悉分形背后的数学知识时开始展开一些自我探索，同时我们将使用功能强大的 NumPy 通用函数来高效执行必要的计算。

您的操作

• 编写函数以绘制各种 Julia 集

• 创建 Mandelbrot 集的可视化效果

• 编写一个函数以计算 Newton 分形

• 尝试各种通用分形类型

您将了解

• 更好地直观了解分形在数学上的工作原理

• 对 NumPy 通用函数和布尔索引的基本了解

• 在 NumPy 中处理复数的基础知识

• 如何创建您自己的独特分形可视化效果

您需要

• Matplotlib

• make_axis_locatable 函数来自 mpl_toolkits API

可以如下导入

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

• 熟悉一些 Python、NumPy 和 matplotlib

• 了解初等数学函数，例如 指数、 正弦、多项式 等

• 非常基本地了解 复数 将很有用

• 了解 导数 可能会很有帮助

热身

为了直观地了解分形是什么，我们将从一个示例开始。

考虑以下等式

\(f(z) = z^2 -1 \)

其中 z 是一个复数（即形式 \(a + bi\) 的复数）

为了方便，我们将为它编写一个 Python 函数

def f(z):
    return np.square(z) - 1

请注意，我们使用的平方函数是 NumPy 通用函数 的一个示例；我们很快将回到此决策的重要性。

为了直观地了解函数的行为，我们可以尝试插入一些不同的值。

对于 \(z = 0\)，我们希望得到 \(-1\)

f(0)

np.int64(-1)

因为我们在设计中使用了通用函数，所以我们可以同时计算多个输入

z = [4, 1-0.2j, 1.6]
f(z)

array([15.  +0.j , -0.04-0.4j,  1.56+0.j ])

一些值增长，一些值收缩，另一些值变化不大。

为了在大范围内查看函数的行为，我们可以将函数应用于复平面的子集并绘制结果。为了创建子集（或网格），我们可以使用 meshgrid 函数。

x, y = np.meshgrid(np.linspace(-10, 10, 20), np.linspace(-10, 10, 20))
mesh = x + (1j * y)  # Make mesh of complex plane

现在我们将会向网格中包含的每个值应用我们的函数。因为我们在设计中使用了一个通用函数，这意味着我们可以一次传入整个网格。出于两个原因，这样极其方便：它减少了所需要编写代码的量，并且极大提高了效率（因为通用函数在其计算中使用了系统级 C 编程）。

在这里，我们在使用 3D 散点图 的函数“迭代”一次后，绘制网格中每个元素的绝对值（或模数）。

output = np.abs(f(mesh))  # Take the absolute value of the output (for plotting)

fig = plt.figure()
ax = plt.axes(projection='3d')

ax.scatter(x, y, output, alpha=0.2)

ax.set_xlabel('Real axis')
ax.set_ylabel('Imaginary axis')
ax.set_zlabel('Absolute value')
ax.set_title('One Iteration: $ f(z) = z^2 - 1$');

这为我们提供了一个关于函数一次迭代所做工作的粗略想法。某些区域（尤其是在最靠近 \((0,0i)\) 的区域）仍然很小，而其他区域增长很大。请注意，我们通过取绝对值失去了有关输出的信息，但这是我们能够绘制图形的唯一方式。

让我们看看当我们向网格应用 2 次迭代时会发生什么。

output = np.abs(f(f(mesh)))

ax = plt.axes(projection='3d')

ax.scatter(x, y, output, alpha=0.2)

ax.set_xlabel('Real axis')
ax.set_ylabel('Imaginary axis')
ax.set_zlabel('Absolute value')
ax.set_title('Two Iterations: $ f(z) = z^2 - 1$');

再一次，我们看到原点附近的仍然很小，而绝对值（或模数）较大的值“爆炸”了。

从第一印象来看，它的行为似乎正常，甚至可能显得平庸。分形往往有超出我们肉眼所见的内容；当我们开始应用更多迭代时，奇异的行为就会显现出来。

考虑三个复数：

\(z_1 = 0.4 + 0.4i \),

\(z_2 = z_1 + 0.1\),

\(z_3 = z_1 + 0.1i\)

鉴于我们前两个图形的形状，我们预计当我们向它们应用迭代时，这些值将保持在原点附近。我们来看看当我们向每个值应用 10 次迭代时会发生什么。

selected_values = np.array([0.4 + 0.4j, 0.41 + 0.4j, 0.4 + 0.41j])
num_iter = 9

outputs = np.zeros((num_iter+1, selected_values.shape[0]), dtype=complex)
outputs[0] = selected_values

for i in range(num_iter):
    outputs[i+1] = f(outputs[i])  # Apply 10 iterations, save each output

fig, axes = plt.subplots(1, selected_values.shape[0], figsize=(16, 6))
axes[1].set_xlabel('Real axis')
axes[0].set_ylabel('Imaginary axis')

for ax, data in zip(axes, outputs.T):
    cycle = ax.scatter(data.real, data.imag, c=range(data.shape[0]), alpha=0.6)
    ax.set_title(f'Mapping of iterations on {data[0]}')

fig.colorbar(cycle, ax=axes, location="bottom", label='Iteration');

令我们惊讶的是，函数的行为与我们的假设完全不符。这是分形具有的混沌行为的一个典型例子。在前两个图形中，该值在最后一次迭代中“爆炸”，远远跳出了它先前所包含的区域。而第三个图形则保持在靠近原点的一个小区域内，尽管值发生了细微的变化，却产生了完全不同的行为。

这就引出了我们一个极其重要的问题：在值发散（“爆炸”）前，对每个值可以应用多少次迭代？

正如我们在前两个图形中看到的那样，值离原点越远，它们爆炸的速度一般就越快。虽然对于较小的值（如 \(z_1, z_2, z_3\)），其行为是不确定的，但我们可以假设，如果一个值超过了从原点到某个距离（比如 2），那么它必然会发散。我们将这个阈值称为半径。

这让我们无需执行过多的计算就能量化给定值下这个函数的行为。一旦超过半径，我们就可以停止迭代，这给了我们一种方法来回答我们提出的问题。如果我们在发散之前统计出了应用了多少计算，我们将深入了解函数的行为，否则很难实现追踪。

当然，我们能够做的更好，设计一个在整个网格上执行过程的函数。

def divergence_rate(mesh, num_iter=10, radius=2):

    z = mesh.copy()
    diverge_len = np.zeros(mesh.shape)  # Keep tally of the number of iterations

    # Iterate on element if and only if |element| < radius (Otherwise assume divergence)
    for i in range(num_iter):
        conv_mask = np.abs(z) < radius
        diverge_len[conv_mask] += 1
        z[conv_mask] = f(z[conv_mask])

    return diverge_len

乍一看，这个函数的行为可能让人困惑，所以，解释一些符号是有帮助的。

我们的目标是对网格中的每个值进行迭代，并在值发散前统计迭代次数。由于一些值会比其他值发散得快，我们需要一个过程，它仅对绝对值足够小的值进行迭代。我们还希望在值一超过半径后，就停止对值进行统计。为此，我们可以使用布尔索引，这是一项 NumPy 特性，与通用函数配对时无可匹敌。布尔索引允许在 NumPy 数组中有条件地执行操作，而无需对每个数组值单独进行循环和检查。

具体来说，我们使用循环将迭代应用于我们的函数 \(f(z) = z^2 -1 \) 并进行统计。使用布尔索引，我们仅对绝对值小于 2 的值应用迭代。

完成这项工作后，我们就可以绘制我们的第一个分形！我们将使用imshow函数来创建统计的彩色编码可视化。

x, y = np.meshgrid(np.linspace(-2, 2, 400), np.linspace(-2, 2, 400))
mesh = x + (1j * y)

output = divergence_rate(mesh)

fig = plt.figure(figsize=(5, 5))
ax = plt.axes()

ax.set_title('$f(z) = z^2 -1$')
ax.set_xlabel('Real axis')
ax.set_ylabel('Imaginary axis')

im = ax.imshow(output, extent=[-2, 2, -2, 2])
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
plt.colorbar(im, cax=cax, label='Number of iterations');

这个惊人的视觉表现传达的是函数行为的复杂性。黄色区域代表保持较小值的区域，而紫色区域代表发散值。出现在收敛值和发散值边界上的美丽图案在你意识到它是从如此简单的函数创建时，会更令人着迷。

朱莉娅集

我们刚才探讨的是特定朱莉娅集的分形可视化的示例。

考虑当 \(c\) 表示复数时函数 \(f(z) = z^2 + c\)。 \(c\) 的填充的朱利亚集是所有复数 z 的集合，其中函数在 \(f(z)\) 收敛。类似地，填充的朱利亚集的边界就是我们称之为朱利亚集的内容。在我们上面的可视化中，我们可以看到黄色区域表示 \(c = -1\) 填充的朱利亚集的近似值，而黄绿色边框将包含朱利亚集。

为了获得更广泛的“朱利亚分形”，我们可以编写一个允许传入 \(c\) 的不同值的函数

def julia(mesh, c=-1, num_iter=10, radius=2):

    z = mesh.copy()
    diverge_len = np.zeros(z.shape)

    for i in range(num_iter):
        conv_mask = np.abs(z) < radius
        z[conv_mask] = np.square(z[conv_mask]) + c
        diverge_len[conv_mask] += 1

    return diverge_len

为了使我们的生活更轻松，我们将创建几个网格，这些网格将在其余示例中重复使用

x, y = np.meshgrid(np.linspace(-1, 1, 400), np.linspace(-1, 1, 400))
small_mesh = x + (1j * y)

x, y = np.meshgrid(np.linspace(-2, 2, 400), np.linspace(-2, 2, 400))
mesh = x + (1j * y)

我们还将编写一个函数，该函数将用于创建我们的分形图

def plot_fractal(fractal, title='Fractal', figsize=(6, 6), cmap='rainbow', extent=[-2, 2, -2, 2]):

    plt.figure(figsize=figsize)
    ax = plt.axes()

    ax.set_title(f'${title}$')
    ax.set_xlabel('Real axis')
    ax.set_ylabel('Imaginary axis')

    im = ax.imshow(fractal, extent=extent, cmap=cmap)
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.1)
    plt.colorbar(im, cax=cax, label='Number of iterations')

使用我们新定义的函数，我们可以再次快速绘制第一个分形

output = julia(mesh, num_iter=15)
kwargs = {'title': 'f(z) = z^2 -1'}

plot_fractal(output, **kwargs);

我们还可以通过尝试 \(c\) 的不同值来探索一些不同的朱利亚集。令人惊讶的是它对分形的形状有多大影响。

例如，设置 \(c = \frac{\pi}{10}\) 会给我们一个非常优雅的云朵形状，而设置 c = \(-\frac{3}{4} + 0.4i\) 会产生一个完全不同的图案。

output = julia(mesh, c=np.pi/10, num_iter=20)
kwargs = {'title': r'f(z) = z^2 + \dfrac{\pi}{10}', 'cmap': 'plasma'}

plot_fractal(output, **kwargs);

output = julia(mesh, c=-0.75 + 0.4j, num_iter=20)
kwargs = {'title': r'f(z) = z^2 - \dfrac{3}{4} + 0.4i', 'cmap': 'Greens_r'}

plot_fractal(output, **kwargs);

曼德布罗特集

与朱利亚集密切相关的是著名的曼德布罗特集，其定义略有不同。再次，我们定义 \(f(z) = z^2 + c\)，其中 \(c\) 是一个复数，但这一次我们的重点是选择 \(c\)。我们说，如果 f 在 \(z = 0\) 收敛，那么 \(c\) 是曼德布罗特集的元素。一个等效的定义是，如果对 \(f(c)\) 进行无限次迭代并且不会“爆炸”，那么 \(c\) 是曼德布罗特集的元素。我们将稍稍调整我们的朱利亚函数（并相应地重命名它），以便我们可以绘制具有优雅分形图案的曼德布罗特集的可视化效果。

def mandelbrot(mesh, num_iter=10, radius=2):

    c = mesh.copy()
    z = np.zeros(mesh.shape, dtype=np.complex128)
    diverge_len = np.zeros(z.shape)

    for i in range(num_iter):
        conv_mask = np.abs(z) < radius
        z[conv_mask] = np.square(z[conv_mask]) + c[conv_mask]
        diverge_len[conv_mask] += 1

    return diverge_len

output = mandelbrot(mesh, num_iter=50)
kwargs = {'title': 'Mandelbrot \\ set', 'cmap': 'hot'}

plot_fractal(output, **kwargs);

推广朱利亚集

通过为通用的函数提供参数，我们可以更进一步地推广我们的 Julia 函数，让我们可以传入该函数。这将允许我们绘制形式为 \(f(z) = g(z) + c\) 的分形，其中 g 是我们选定的通用函数。

def general_julia(mesh, c=-1, f=np.square, num_iter=100, radius=2):

    z = mesh.copy()
    diverge_len = np.zeros(z.shape)

    for i in range(num_iter):
        conv_mask = np.abs(z) < radius
        z[conv_mask] = f(z[conv_mask]) + c
        diverge_len[conv_mask] += 1

    return diverge_len

可以使用我们的通用 Julia 函数绘制的一组很酷的分形是形式为 \(f(z) = z^n + c\) 的分形，其中 \(n\) 是正整数。出现了一个非常酷的模式，即“突出”区域的数量与我们在迭代时将函数提高的次数相匹配。

fig, axes = plt.subplots(2, 3, figsize=(8, 8))
base_degree = 2

for deg, ax in enumerate(axes.ravel()):
    degree = base_degree + deg
    power = lambda z: np.power(z, degree)  # Create power function for current degree

    diverge_len = general_julia(mesh, f=power, num_iter=15)
    ax.imshow(diverge_len, extent=[-2, 2, -2, 2], cmap='binary')
    ax.set_title(f'$f(z) = z^{degree} -1$')

不用说，可以通过摆弄输入函数、\(c\) 的值、迭代次数、半径甚至网格的密度和颜色的选择来进行大量的探索。

牛顿分形

牛顿分形是一类特定分形，其中迭代包括对函数（通常是多项式）及其导数的比值进行加法或减法，然后将其添加到输入值中。在数学上，它可以表示为

\(z := z - \frac{f(z)}{f'(z)}\)

我们将定义分形的通用版本，该版本允许通过传入我们选择的函数来绘制不同的变化。

def newton_fractal(mesh, f, df, num_iter=10, r=2):

    z = mesh.copy()
    diverge_len = np.zeros(z.shape)

    for i in range(num_iter):
        conv_mask = np.abs(z) < r
        pz = f(z[conv_mask])
        dp = df(z[conv_mask])
        z[conv_mask] = z[conv_mask] - pz/dp
        diverge_len[conv_mask] += 1

    return diverge_len

现在我们可以尝试使用一些不同的函数了。对于多项式，我们可以使用 NumPy 多项式类轻松创建我们的绘图，该类具有用于计算导数的内置功能。

例如，我们尝试一个高次多项式

p = np.polynomial.Polynomial([-16, 0, 0, 0, 15, 0, 0, 0, 1])
p

\[x \mapsto \text{-16.0}\color{LightGray}{ + \text{0.0}\,x}\color{LightGray}{ + \text{0.0}\,x^{2}}\color{LightGray}{ + \text{0.0}\,x^{3}} + \text{15.0}\,x^{4}\color{LightGray}{ + \text{0.0}\,x^{5}}\color{LightGray}{ + \text{0.0}\,x^{6}}\color{LightGray}{ + \text{0.0}\,x^{7}} + \text{1.0}\,x^{8}\]

其导数为

p.deriv()

\[x \mapsto \color{LightGray}{\text{0.0}}\color{LightGray}{ + \text{0.0}\,x}\color{LightGray}{ + \text{0.0}\,x^{2}} + \text{60.0}\,x^{3}\color{LightGray}{ + \text{0.0}\,x^{4}}\color{LightGray}{ + \text{0.0}\,x^{5}}\color{LightGray}{ + \text{0.0}\,x^{6}} + \text{8.0}\,x^{7}\]

output = newton_fractal(mesh, p, p.deriv(), num_iter=15, r=2)
kwargs = {'title': r'f(z) = z - \dfrac{(z^8 + 15z^4 - 16)}{(8z^7 + 60z^3)}', 'cmap': 'copper'}

plot_fractal(output, **kwargs)

漂亮！让我们再尝试一个

f(z) = \(tan^2(z)\)

\(\frac{df}{dz} = 2 \cdot tan(z) sec^2(z) =\frac{2 \cdot tan(z)}{cos^2(z)}\)

这使得 \(\frac{f(z)}{f'(z)} = tan^2(z) \cdot \frac{cos^2(z)}{2 \cdot tan(z)} = \frac{tan(z)\cdot cos^2(z)}{2} = \frac{sin(z)\cdot cos(z)}{2}\)

def f_tan(z):
    return np.square(np.tan(z))


def d_tan(z):
    return 2*np.tan(z) / np.square(np.cos(z))

output = newton_fractal(mesh, f_tan, d_tan, num_iter=15, r=50)
kwargs = {'title': r'f(z) = z - \dfrac{sin(z)cos(z)}{2}', 'cmap': 'binary'}

plot_fractal(output, **kwargs);

请注意，有时你必须调整半径才能得到一个外观不错的分形。

最后，我们可以对选择函数时放肆一点

\(f(z) = \sum_{i=1}^{10} sin^i(z)\)

\(\frac{df}{dz} = \sum_{i=1}^{10} i \cdot sin^{i-1}(z) \cdot cos(z)\)

def sin_sum(z, n=10):
    total = np.zeros(z.size, dtype=z.dtype)
    for i in range(1, n+1):
        total += np.power(np.sin(z), i)
    return total


def d_sin_sum(z, n=10):
    total = np.zeros(z.size, dtype=z.dtype)
    for i in range(1, n+1):
        total += i * np.power(np.sin(z), i-1) * np.cos(z)
    return total

我们将称这个为“出格分形”，因为它的方程不好看，无法放在标题中。

output = newton_fractal(small_mesh, sin_sum, d_sin_sum, num_iter=10, r=1)
kwargs = {'title': 'Wacky \\ fractal', 'figsize': (6, 6), 'extent': [-1, 1, -1, 1], 'cmap': 'terrain'}

plot_fractal(output, **kwargs)

令人惊叹的是，这些分形彼此既不同又相似。这让我们步入最后一个部分。

创建自己的分形

一旦你熟悉了基础知识，分形就会变得更加令人兴奋，因为还有很多内容可以探索。现在，我们将通过探索创建唯一分形的一些不同方式来结束我们的教程。我鼓励你自己尝试一些东西（如果你还没有这样做的话）。

可以尝试的第一个地方是广义朱丽亚集合的函数，我们可以尝试将不同的函数作为参数传递进去。

让我们从选择开始

\(f(z) = tan(z^2)\)

def f(z):
    return np.tan(np.square(z))

output = general_julia(mesh, f=f, num_iter=15, radius=2.1)
kwargs = {'title': 'f(z) = tan(z^2)', 'cmap': 'gist_stern'}

plot_fractal(output, **kwargs);

如果我们将定义的函数组合在一个正弦函数中，会发生什么？

让我们尝试定义

\(g(z) = sin(f(z)) = sin(tan(z^2))\)

def g(z):
    return np.sin(f(z))

output = general_julia(mesh, f=g, num_iter=15, radius=2.1)
kwargs = {'title': 'g(z) = sin(tan(z^2))', 'cmap': 'plasma_r'}

plot_fractal(output, **kwargs);

接下来，让我们创建一个函数，在每次迭代中将 f 和 g 应用于输入并相加结果

\(h(z) = f(z) + g(z) = tan(z^2) + sin(tan(z^2))\)

def h(z):
    return f(z) + g(z)

output = general_julia(small_mesh, f=h, num_iter=10, radius=2.1)
kwargs = {'title': 'h(z) = tan(z^2) + sin(tan(z^2))', 'figsize': (7, 7), 'extent': [-1, 1, -1, 1], 'cmap': 'jet'}

plot_fractal(output, **kwargs);

你甚至可以通过自己的错误创建漂亮的图形。这里有一个错误地计算牛顿分形导数而意外创建的图形

def accident(z):
    return z - (2 * np.power(np.tan(z), 2) / (np.sin(z) * np.cos(z)))

output = general_julia(mesh, f=accident, num_iter=15, c=0, radius=np.pi)
kwargs = {'title': 'Accidental \\ fractal', 'cmap': 'Blues'}

plot_fractal(output, **kwargs);

不用说，只要摆弄各种 NumPy 通用函数的组合并调整参数，就可以创造出无穷无尽的有趣分形作品。

总之

我们今天学习了很多关于生成分形的事情。我们看到如何使用通用函数有效地计算需要多次迭代的复杂分形。我们还利用了布尔索引，这样就可以减少计算量，而无需单独验证每个值。最后，我们了解了很多分形本身。作为回顾：

• 分形图像通过对一组值进行函数迭代创建，并记录每个值需要多长时间才能通过某个阈值

• 图像中的颜色对应于值的计数

• 填充的朱莉亚集合 \(c\) 由其中 \(f(z) = z^2 + c\) 收敛的所有复数 z 组成

• 朱莉亚集合 \(c\) 是构成填充朱莉亚集合边界的复数的集合

• 曼德尔布罗特集合是其中 \(f(z) = z^2 + c\) 在 0 收敛的所有值 \(c\)

• 牛顿分形使用以下形式的函数：\(f(z) = z - \frac{p(z)}{p'(z)}\)

• 当你调整迭代次数、收敛半径、网格尺寸、颜色、函数选择和参数选择时，分形图像可能会有所不同

自行操作

• 使用广义朱丽娅集函数的参数进行尝试，尝试使用常数值、迭代次数、函数选择、半径和颜色选择进行尝试。

• 访问维基百科上的“按豪斯多夫维数列出的分形列表”（在延伸阅读部分中提供链接）页面，尝试编写本教程中未提及的分形的函数。

延伸阅读

• 有关分形理论的更多信息

• 有关朱丽娅集的延伸阅读

• 有关曼德尔布罗集的更多详细信息

• 牛顿分形的更完整处理

• 各种分形的列表