绘制分形图#
分形是一种美丽而引人入胜的数学形式，通常可以通过相对简单的指令集来创建。在自然界中，它们存在于各种地方，如海岸线、贝壳和蕨类植物，甚至被用于制造某些类型的天线。分形的数学思想早已为人所知，但直到20世纪70年代，随着计算机图形学的进步和一些意外的发现，研究人员如Benoît Mandelbrot才偶然发现了分形所拥有的真正神秘的视觉效果，它们才开始受到真正的重视。
今天，我们将学习如何绘制这些美丽的视觉效果图，并随着我们对分形背后的数学原理越来越熟悉，开始进行一些探索。我们将使用强大的NumPy通用函数（universal functions）来高效地执行必要的计算。
你将做什么#
编写一个函数来绘制各种Julia集
创建Mandelbrot集的视觉效果图
编写一个计算牛顿分形的函数
实验通用分形类型的变体
你将学到什么#
更好地直观理解分形在数学上的运作方式
对NumPy通用函数和布尔索引（Boolean Indexing）的基本理解
在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通用函数（universal function）的一个示例；我们稍后将回到这一决定的重要性。
为了直观理解该函数的行为，我们可以尝试代入一些不同的值。
对于\(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
这个函数的行为乍一看可能令人困惑，所以解释一些符号会有所帮助。
我们的目标是迭代网格中的每个值，并统计值发散之前的迭代次数。由于某些值会比其他值更快地发散，我们需要一个只迭代绝对值足够小的值的过程。我们还希望一旦值超过半径就停止计数。为此，我们可以使用布尔索引（Boolean Indexing），这是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集（filled-in Julia set）是函数在\(f(z)\)处收敛的所有复数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\)。如果函数在\(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多项式类（Polynomial class）来创建图形，它内置了计算导数的功能。
例如，让我们尝试一个更高次的多项式
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通用函数的各种组合以及调整参数，就可以创造出几乎无限多的有趣分形。
总结#
今天我们学到了很多关于生成分形的知识。我们看到了如何使用通用函数高效地计算需要多次迭代的复杂分形。我们还利用了布尔索引，这使得在无需单独验证每个值的情况下减少了计算量。最后，我们对分形本身也有了深入了解。回顾一下
分形图像是通过对一组值迭代一个函数，并记录每个值通过特定阈值所需的时间来创建的
图像中的颜色对应于值的计数
对于\(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维度分类的分形列表”（在“延伸阅读”部分链接），并尝试为本教程中未提及的分形编写一个函数。