绘制分形#
分形是美丽的、引人入胜的数学形式，通常可以通过相对简单的指令集创建。在自然界中，它们可以在各种地方找到，如海岸线、海螺和蕨类植物，甚至被用于创建某些类型的天线。分形的数学概念早已为人所知，但随着计算机图形学的进步和一些偶然的发现，研究人员如Benoît Mandelbrot在1970年代开始真正欣赏到分形所拥有的真正神秘的视觉效果。
今天我们将学习如何绘制这些美丽的可视化效果，并开始为自己做一些探索，因为我们熟悉了分形背后的数学，并将使用功能强大的 NumPy 通用函数来高效地执行必要的计算。
你要做什么#
编写一个用于绘制各种 Julia 集的函数
创建一个Mandelbrot集的可视化
编写一个计算牛顿分形的函数
实验不同类型的广义分形变体
你将学到什么#
更好地理解分形在数学上是如何工作的
关于NumPy通用函数和布尔索引的基本理解
在 NumPy 中使用复数的基础
如何创建你自己独特的分形可视化
你需要什么#
mpl_toolkits API 中的
make_axis_locatable函数
可以按如下方式导入：
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
热身#
为了对分形是什么有一些直觉，我们将从一个例子开始。
考虑以下方程：
\(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\) 是一个复数。填充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
它有导数：
p.deriv()
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
我们将这个称为 ‘Wacky fractal’，因为它的方程放在标题里会不太有趣。
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 通用函数组合和调整参数，可以创造出几乎无穷无尽的有趣分形作品。
总之#
今天我们学到了很多关于生成分形的内容。我们看到了如何使用通用函数高效地计算需要多次迭代才能生成的复杂分形。我们还利用了布尔索引，这使得在不逐个验证每个值的情况下减少了计算量。最后，我们学到了很多关于分形本身的知识。总结一下：
分形图像通过在一组值上迭代一个函数来创建，并记录每个值超过某个阈值所需的时间。
图像中的颜色对应于值的计数
填充的 Julia 集对于 \(c\) 包含所有复数
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 集函数的参数，尝试调整常数值、迭代次数、函数选择、半径和颜色选择。
访问“分形列表按豪斯多夫维数”维基百科页面（在进一步阅读部分链接）并尝试为一个未在此教程中提到的分形编写一个函数。