MENU

Numpy でマンデルブロ集合を作成する方法について

NumPy
NumPy Python
  1. ホーム
  2. Python
  3. NumPy
  4. Numpy でマンデルブロ集合を作成する方法について
目次
  1. 概要
  2. マンデルブロ集合
  3. マンデルブロ集合の性質
  4. Numpy による実装例
    1. マンデルブロ集合を作成する
    2. マンデルブロ集合を白黒画像で可視化する
    3. マンデルブロ集合をカラー画像で可視化する
    4. マンデルブロ集合の一部を拡大する

概要

Python の数値計算ライブラリ NumPy を使って、フラクタルであるマンデルブロ集合 (Mandelbrot set) を作成し、画像で可視化する方法を解説します。

マンデルブロ集合

{z0(c)=0zn+1 (c)=(zn(c))2+c \begin{cases} z_0(c) = 0 \\ z_{n + 1}\ (c) = (z_n(c))^2 + c \end{cases}

で定義される複素数列 {zn(c)}\{z_n(c)\} が n→∞n \to \infty で無限大に発散しない条件を満たす複素数 c∈Cc \in \mathbb{C} 全体の集合をマンデルブロ集合 (Mandelbrot set) といい、M\mathbb{M} で表します。

例:

  • c=1c = 1 のとき、数列 {zn}\{z_n\} は 0,1,2,5,26,⋯0, 1, 2, 5, 26, \cdots となり発散するので、マンデルブロ集合に含まれます。
  • c=−1c = -1 のとき、数列 {zn}\{z_n\} は 0,−1,0,−1,0,⋯0, -1, 0, -1, 0, \cdots となり発散しないので、マンデルブロ集合に含まれません。

マンデルブロ集合の性質

コンピューターでは、極限 lim⁡n→∞zn\displaystyle \lim_{n \to \infty} z_n は計算できないため、発散するかどうかをどのように判定するかが問題になります。この問題に以下の定理が活用できます。

c∈Cc\in \mathbb{C} に対して、∣zn(c)∣>2|z_n(c)| > 2 となる自然数 nn が存在すれば、複素数列 {zn(c)}\{z_n(c)\} は発散する。すなわち、cc はマンデルブロ集合に含まれない。 逆に cc がマンデルブロ集合に含まれなければ、∣zn(c)∣>2|z_n(c)| > 2 となる自然数 nn が存在する。

証明は以下を参照してください。

problem solving – Show that Mandelbrot set is contained within the closed disc of r=2 – Mathematics Stack Exchange

この定理により、第 nn 項まで計算した時点で ∣zn(c)∣>2|z_n(c)| > 2 となれば、cc はマンデルブロ集合に含まれないと判定できます。

Numpy による実装例

マンデルブロ集合を作成する

先にマンデルブロ集合を作成するコードを示します。

In [1]: 
import numpy as np
from PIL import Image


def generate_mandelbrot(
    real_range=(-2.0, 2.0),
    imag_range=(-2.0, 2.0),
    resolution=(400, 400),
    iterations=100,
):
    """マンデルブロ集合を作成する。

    Keyword Arguments:
        real_range: 実軸の範囲を (下限, 上限) のタプルで指定する。
        imag_range: 虚軸の範囲を (下限, 上限) のタプルで指定する。
        resolution: 解像度を (実軸の解像度, 虚軸の解像度) のタプルで指定する。
        max_iters: 反復回数
    """
    # 複素平面上に格子状の点を作成する。
    R, I = np.meshgrid(
        np.linspace(*real_range, resolution[0]),
        np.linspace(*imag_range, resolution[1]),
    )
    C = R + I * 1j

    # 複素数列 z_i(c)
    Z = np.zeros_like(C)
    # 発散が確定するまでの反復回数 (|z_i(c)| <= 2 を満たす最大の i)
    step = np.zeros(C.shape, dtype=int)

    for i in range(iterations):
        # 複素数列 z_i(c) を計算する。|z_i(c)| > 2 の場合、オーバーフローを防ぐために更新しない。
        Z = np.where(np.abs(Z) <= 2, Z ** 2 + C, Z)

        # |z_i(c)| <= 2 なら +1 する。
        step += np.abs(Z) <= 2

    return Z, step


# マンデルブロ集合を作成する。
real_range = (-2.0, 2.0)  # 実軸の範囲
imag_range = (-2.0, 2.0)  # 虚軸の範囲
resolution = (500, 500)  # (実軸の解像度, 虚軸の解像度)
iterations = 300  # 反復回数
Z, step = generate_mandelbrot(real_range, imag_range, resolution, iterations)
Python

実軸、虚軸の範囲をそれぞれ real_range、imag_range で指定します。 複素平面のこの範囲に解像度 resolution で格子状の点を作成します。

real_range = (-2.0, 2.0)
imag_range = (-2.0, 2.0)
resolution = (10, 10)
R, I = np.meshgrid(
    np.linspace(*real_range, resolution[0]), np.linspace(*imag_range, resolution[1]),
)
Python

for 内では、複素数列 zi(c)z_i(c) を計算します。 zi+1(c)=(zi(c))2+cz_{i + 1}(c) = (z_i(c))^2 + c で数列を更新していきますが、∣zi(c)∣>2|z_i(c)| > 2 の場合は発散することは確定しているので、更新を続けた結果、オーバーフローしてしまうのを防ぐために更新を止めます。

また、∣zi(c)∣≤2|z_i(c)| \le 2 を満たす最大の ii も計算します。 この計算には、True は1、False は0と同値であることを利用し、比較 ∣zi(c)∣≤2|z_i(c)| \le 2 の結果である bool 値を足していきます。 この情報はマンデルブロ集合を色付きで画像化する際に利用します。

for i in range(iterations):
    # 複素数列 z_i(c) を計算する。|z_i(c)| > 2 の場合、オーバーフローを防ぐために更新しない。
    Z = np.where(np.abs(Z) <= 2, Z ** 2 + C, Z)

    # |z_i(c)| <= 2 なら +1 する。
    step += np.abs(Z) <= 2
Python

マンデルブロ集合を白黒画像で可視化する

マンデルブロ集合を白黒画像で画像化する例を示します。

∣zi(c)∣≤2|z_i(c)| \le 2 である複素数 cc は、発散していないので黒色 (0)、そうでないなら発散しているので白色 (255) で画像にします。

In [2]: 
def simple_colorize(Z):
    # マンデルブロ集合は黒、それ以外は白として2値画像を作成する。
    img = np.where(np.abs(Z) <= 2, 0, 255).astype(np.uint8)
    img = Image.fromarray(img)

    return img


simple_colorize(Z)
Python

マンデルブロ集合をカラー画像で可視化する

マンデルブロ集合をカラー画像で画像化する例を示します。

マンデルブロ集合に含まれる点 (∣zi(c)∣≤2|z_i(c)| \le 2 である点 cc) の色は、前項と同様に黒色 (0,0,0)(0, 0, 0) にします。 発散している点も発散が確定した反復回数によって色を変えます。step / iterations で [0,1][0, 1] の範囲の値にし、matplotlib のカラーマップで色に変換します。

In [3]: 
def colorize(Z, step, iterations, cmap_name="jet"):
    # 各点を発散が確定した反復回数に応じて色分けする。
    cmap = plt.cm.get_cmap(cmap_name)
    img = cmap(step / iterations, bytes=True)

    # マンデルブロ集合に含まれる点の色は黒にする。
    img[np.abs(Z) <= 2] = 0
    img = Image.fromarray(img)

    return img


img = colorize(Z, step, iterations)
img
Python

cmap_name に matplotlib のカラーマップ名を指定し、カラーマップを変更すると、また違った雰囲気の画像が作成できます。

マンデルブロ集合の一部を拡大する

マンデルブロ集合の一部を拡大して描画したい場合は、real_range、imag_range の範囲を狭くします。 以下、c=−0.75+0.1ic = -0.75 + 0.1i の周辺を画像化した例になります。 タツノオトシゴのような模様があるので、sea horse valley (タツノオトシゴの谷) と呼ばれています。

In [4]: 
# マンデルブロ集合の -0.75 + 0.1i 近辺の図
center = -0.75 + 0.1j
offset = 0.005

real_range = (center.real - offset, center.real + offset)
imag_range = (center.imag - offset, center.imag + offset)
resolution = (500, 500)
iterations = 1000
Z, step = generate_mandelbrot(real_range, imag_range, resolution, iterations)

img = colorize(Z, step, iterations)
img
Python
NumPy
NumPy Python

コメント

コメントする

目次
  1. 概要
  2. マンデルブロ集合
  3. マンデルブロ集合の性質
  4. Numpy による実装例
    1. マンデルブロ集合を作成する
    2. マンデルブロ集合を白黒画像で可視化する
    3. マンデルブロ集合をカラー画像で可視化する
    4. マンデルブロ集合の一部を拡大する