绘制分形图

分形是美丽而引人入胜的数学形式，通常可以从一组相对简单的指令中创建出来。在自然界中，它们存在于各种地方，例如海岸线、贝壳和蕨类植物，甚至被用于创建某些类型的触角。分形的数学思想已经存在相当一段时间了，但直到 20 世纪 70 年代，随着计算机图形学的进步和一些偶然的发现，像Benoît Mandelbrot这样的研究人员才真正开始欣赏分形所拥有的真正令人费解的可视化效果。

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

你将做什么

• 编写一个用于绘制各种 Julia 集的函数

• 创建 Mandelbrot 集的可视化效果

• 编写一个计算牛顿分形的函数

• 尝试各种一般分形类型的变化

你将学到什么

• 对分形如何从数学上运作有更好的直觉

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

• 在 NumPy 中使用复数的基础知识

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

你需要什么

• Matplotlib

• 来自 mpl_toolkits API 的make_axis_locatable 函数

可以按如下方式导入

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');

这张令人惊叹的视觉图展现了函数行为的复杂性。黄色区域代表保持较小的值，而紫色区域代表发散的值。收敛值和发散值边界上出现的美丽图案，当你意识到它是从如此简单的函数中产生的时，就更加令人着迷了。

Julia集

我们刚才探索的是特定Julia集的分形可视化示例。

考虑函数\(f(z) = z^2 + c\)，其中\(c\)是复数。 \(c\)的填充Julia集是所有复数z的集合，其中函数在\(f(z)\)处收敛。同样，填充Julia集的边界就是我们所说的Julia集。在我们上面的可视化中，我们可以看到黄色区域代表\(c = -1\)的填充Julia集的近似值，而黄绿色的边界包含Julia集。

为了访问更广泛的“Julia分形”，我们可以编写一个函数，允许传入不同的\(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\)值来探索一些不同的Julia集。它对分形的形状的影响之大令人惊讶。

例如，设置\(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);

Mandelbrot集

与Julia集密切相关的是著名的Mandelbrot集，它的定义略有不同。我们再次定义\(f(z) = z^2 + c\)，其中\(c\)是复数，但这次我们的重点是我们对\(c\)的选择。如果f在\(z = 0\)处收敛，我们就说\(c\)是Mandelbrot集的元素。等效的定义是说，如果\(f(c)\)可以无限迭代而不“爆炸”，则\(c\)是Mandelbrot集的元素。我们将稍微调整我们的Julia函数（并将其重命名为适当的名称），以便我们可以绘制Mandelbrot集的可视化图像，它具有优雅的分形图案。

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集的推广

我们可以通过为其提供一个参数来进一步推广我们的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)

这些分形彼此之间是如此独特却又如此相似，这真是令人着迷。这将我们引向最后一部分。

创建你自己的分形

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

第一个实验的地方将是广义Julia集的函数，我们可以尝试传入不同的函数作为参数。

让我们从选择开始

\(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\)的填充Julia集包含所有复数z，其中\(f(z) = z^2 + c\)收敛。

• \(c\)的Julia集是构成填充Julia集边界的复数集。

• Mandelbrot集是所有\(c\)值，其中\(f(z) = z^2 + c\)在0处收敛。

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

• 随着你调整迭代次数、收敛半径、网格大小、颜色、函数选择和参数选择，分形图像也会发生变化。

自己动手

• 尝试调整广义Julia集函数的参数，尝试调整常数值、迭代次数、函数选择、半径和颜色选择。

• 访问维基百科的“按Hausdorff维数列出的分形”页面（在“进一步阅读”部分链接），并尝试为本教程中未提及的分形编写函数。

进一步阅读

• 关于分形背后理论的更多信息

• 关于Julia集的进一步阅读

• 关于Mandelbrot集的更多细节

• 对牛顿分形的更完整处理

• 不同分形的列表