跳轉到

8.3 曼德博集合#

概念說明#

曼德博集合（Mandelbrot set）是複平面上最著名的分形，定義只有一句話：

\[ c \in M \iff \text{迭代 } z_{k+1} = z_k^2 + c,\ z_0 = 0 \text{ 的軌道有界} \]

一個平方、一個加法，反覆套用，卻在複平面刻出無限複雜的邊界——放大任何角落都湧現新的自我相似細節，且邊界的分形維數為 2。數學上： \(c \in M\) 當且僅當對應的 Julia 集合連通；曼德博集合是所有 Julia 集合的「目錄」。

計算上無法直接判定「有界」，實用判準是逃逸半徑：若 \(|z_k| > 2\) ，則軌道必然逃至無窮（因 \(|z_{k+1}| \ge |z_k|^2 - |c| > 2|z_k| - |c| \to \infty\) ），可提前停機；若迭代到上限（如 500 步）仍未逃逸，就判為屬於 \(M\) 。這是逃逸時間演算法：對每個像素執行迭代、數逃逸步數、按步數上色。它把「收斂性分析」反轉成繪圖工具——我們不問迭代收斂到哪，只問它多快發散。

工程要點：

  1. 複數迭代極簡： \(z^2 + c\) 用實虛部寫成 \(x \leftarrow x^2 - y^2 + c_x\) ， \(y \leftarrow 2xy + c_y\) ，無三角函數。
  2. 平滑著色：逃逸步數是整數，直接上色會有色帶；用 \(n + 1 - \log_2 \frac{\ln|z_n|}{\ln 2}\) 得連續值。
  3. 效能：每像素最多 \(L\) 次迭代，全圖 \(O(W \cdot H \cdot L)\) ；邊界附近的像素最耗時。加速手法：早期退出、週期檢測、perturbation 理論（深度縮放）。
  4. 停止條件即精度：迭代上限 \(L\) 決定邊界細節的解析度—— \(L\) 太小，集合邊界被「啃掉」；這是「迭代次數即計算資源」最直觀的示範。

曼德博集合的意義超出美學：它是非對稱資訊的縮影——判定「在集合內」沒有有限證書（只能迭代到上限），判定「不在」卻有（逃逸）——一個天然的 NP vs co-NP 直覺教具。

生活類比#

重複「拿利息再投資」：本金小、報酬率適中則財富盤旋在有限範圍（在集合內）；複利一旦過猛就指數爆炸（逃逸）。每個起點 \(c\) 對應一種命運，把所有命運畫在地圖上，邊界正是「穩定與失控」的針鋒相對處——而那裡的形狀無限精細。

演算法步驟#

逃逸時間演算法：

  1. 對視窗內每個像素 \(c = c_x + c_y i\) （ \(x \in [-2.0, 0.8]\) ， \(y \in [-1.4, 1.4]\) ）。
  2. \(z = 0\) ；迭代 \(k = 1 \dots L\) ： \(z \leftarrow z^2 + c\) 。
  3. 若 \(|z| > 2\) ：記錄 \(k\) 為逃逸時間，停。
  4. 若 \(k = L\) 仍未逃逸：判為集合內（通常畫黑色）。
  5. 依 \(k\) （或平滑版）映射到色盤。

Python 範例#

import numpy as np

def mandelbrot(w=400, h=300, max_iter=200,
               x_range=(-2.0, 0.8), y_range=(-1.3, 1.3)):
    xs = np.linspace(*x_range, w)
    ys = np.linspace(*y_range, h)
    C = xs[None, :] + 1j * ys[:, None]     # 複數網格
    Z = np.zeros_like(C)
    escape = np.full(C.shape, max_iter, dtype=int)
    mask = np.ones(C.shape, dtype=bool)    # 尚未逃逸
    for k in range(1, max_iter + 1):
        Z[mask] = Z[mask] * Z[mask] + C[mask]
        newly = mask & (np.abs(Z) > 2.0)
        escape[newly] = k
        mask &= ~newly
        if not mask.any():
            break
    return escape

img = mandelbrot(w=200, h=150, max_iter=100)
inside = np.sum(img == 100)
print(f"內部像素估計: {inside} / {img.size} = {inside/img.size:.1%}")

# 平滑著色值（供渲染）
def smooth(escape, max_iter):
    out = escape.astype(float)
    out[escape < max_iter] += 1 - np.log2(np.log(np.abs(escape[escape < max_iter]) + 1) /
                                          np.log(max_iter))
    return out

# 渲染（需 matplotlib）：
# import matplotlib.pyplot as plt
# plt.imshow(smooth(mandelbrot(800, 600, 500), 500), cmap="magma",
#            extent=(-2.0, 0.8, -1.3, 1.3)); plt.show()

複雜度與收斂性分析#

每像素成本 \(\le L\) 次迭代，全圖 \(O(WH L)\) ，且實際成本集中在邊界帶（內部與遠外部都快速停機）。迭代上限 \(L\) 與可信放大倍率的關係：邊界寬度約 \(O(1/L)\) ，故放大 \(10^4\) 倍需 \(L \sim 10^4\) 。數值誤差隨深度增長，超深縮放需浮點擴充（double-double 或 fixed-point）。相關分析工具：勢函數（green's function） \(G(c) = \lim \frac{\ln|z_k|}{2^k}\) 與距離估計式 \(d(c, M) \approx \frac{|z_k| \ln|z_k|}{|z_k'|}\) （配導數迭代 \(z' \leftarrow 2 z z' + 1\) ），可大幅改善邊界渲染。本節核心訓練：把「會不會收斂」的定性問題，轉化為可執行、可視覺化、可量化的計算。

練習#

  1. 寫出 \(z \leftarrow z^2 + c\) 的實虛部展開，驗證與複數運算結果一致。
  2. 對 \(c = -1\) （週期 2： \(0 \to -1 \to 0\) ）與 \(c = 0.3\) （不動點）手算數步，理解集合內部的動態分類。
  3. 測試 max_iter = 20、100、1000 三種設定，比較邊界的視覺差異（可印出 ASCII 圖），並解釋為何上限不足會「啃掉」邊界。