繪製碎形

碎形是美麗且引人入勝的數學形式，通常可以從相對簡單的一組指令中創建出來。在自然界中，它們可以在許多地方找到，例如海岸線、貝殼和蕨類植物，甚至被用於創建某些類型的天線。碎形的數學概念已經存在相當長一段時間，但直到 1970 年代，隨著電腦圖形技術的進步和一些意外的發現，才真正開始受到重視，引導像 Benoît Mandelbrot 這樣的研究人員偶然發現了碎形所擁有的真正神秘的可視化效果。

今天，我們將學習如何繪製這些美麗的可視化效果，並將開始親自探索，因為我們將熟悉碎形背後的數學原理，並將使用功能強大的 NumPy 通用函數來有效地執行必要的計算。

您將做什麼

• 編寫一個用於繪製各種茱莉亞集的函數

• 創建曼德博集的可視化

• 編寫一個計算牛頓碎形的函數

• 實驗通用碎形類型的變化

您將學到什麼

• 對碎形在數學上如何運作有更好的直覺

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

這個令人驚嘆的視覺效果傳達的是函數行為的複雜性。黃色區域表示保持較小的值，而紫色區域表示發散的值。當您意識到它是從如此簡單的函數創建出來的時候，在收斂值和發散值的邊界上產生的美麗圖案更加令人著迷。

茱莉亞集

我們剛才探索的是特定茱莉亞集的碎形可視化示例。

考慮函數 \(f(z) = z^2 + c\)，其中 \(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\) 的選擇。我們說 \(c\) 是曼德博集的一個元素，如果 f 在 \(z = 0\) 處收斂。一個等效的定義是說 \(c\) 是曼德博集的一個元素，如果 \(f(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);

推廣茱莉亞集

我們可以通過為我們的茱莉亞函數提供一個參數來進一步推廣它，我們想要傳入哪個通用函數。這將允許我們繪製 \(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

可以使用我們的通用茱莉亞函數繪製的一組很酷的碎形是 \(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\) 的填充茱莉亞集由所有複數 z 組成，其中 \(f(z) = z^2 + c\) 收斂

• \(c\) 的茱莉亞集是構成填充茱莉亞集邊界的複數集合

• 曼德博集是所有值 \(c\)，其中 \(f(z) = z^2 + c\) 在 0 處收斂

• 牛頓碎形使用 \(f(z) = z - \frac{p(z)}{p'(z)}\) 形式的函數

• 當您調整迭代次數、收斂半徑、網格大小、顏色、函數選擇和參數選擇時，碎形圖像可能會有所不同

自行練習

• 玩轉通用茱莉亞集函數的參數，嘗試調整常數值、迭代次數、函數選擇、半徑和顏色選擇。

• 訪問維基百科頁面「按郝斯多夫維數排列的碎形列表」（在延伸閱讀部分連結），並嘗試為本教程中未提及的碎形編寫一個函數。

延伸閱讀

• 更多關於碎形背後理論的信息

• 更多關於茱莉亞集的延伸閱讀

• 更多關於曼德博集的詳細信息

• 更完整的牛頓碎形處理

• 不同碎形的列表