logo logo

January 05, 2019 15:49

Pythonでマンデルブロ集合を描画

導入

今回はpythonでのグラフ描画のチュートリアルとして、マンデルブロ集合をプロットしてみます。
これを見たマンデルブロ集合を知っているけど自分で出力はしたことない人には、数値計算やハッキングの楽しさを、
そもそもマンデルブロ集合を知らない人には、簡単な数学とコードの知識だけでこれだけ美しいものを自らの手で描けるということに、少しでも感動していただけたら嬉しいなと思っています。私が高校生の時と指導要領がだいぶ変わっているとは思いますが、
複素数と漸化式をかじったことのある高校生といった数学レベルで内容は理解できると思います。

定義

複数数に馴染みがあるなら、マンデルブロ集合の定義は簡単に理解できます。以下の漸化式
\begin{cases}
z _{n+1} = z _{n}^{2} + c \text{ }(c\in \mathbb{Z})\\
z _{0} = 0
\end{cases}

において、$n \rightarrow \infty$で$z_{n}$が発散しないような複数数$c$の集合がマンデルブロ集合です。
複素平面上にこのような$c$を離散的にプロットしていけば一応、マンデルブロ集合の要素を探してくることができるわけですが、
収束の早い遅いで色付けをすることにより、より美しく描写することができます。

方針

数値計算ですので厳密な収束や発散の判定はできません。そこで適当に大きな数字(実装中のM = 10000000)を定め、これまた適当に大きな上限の$n _{\text{max}}$(実装中のmaxItr)を定めます。$|z _{n}|$がMを超えたら発散判定とします。また超えた場合はその時点での$n$を保存しておき、収束の速さの目安とし、色付けに用いることにします。

実装

まずは判定のためのパラメータを保持し、収束判定計算を行うMandelbrotクラスを実装します。

class Mandelbrot:
    def __init__(self, M, maxItr):
        self.M = M
        self.maxItr = maxItr

    def Calc(self, re, im): # re, im: 2d Array 
        x = 0.0
        y = 0.0

        for i in range(self.maxItr):
            if x**2 + y*2 > self.M:
                return i
            x_ = x**2 - y**2 + re
            y_ = 2*x*y + im
            x, y = x_, y_
        return self.maxItr

次に実際に上記クラスを用いて複素平面上のグリッドに対応する2次元配列に、計算結果を格納していきます。
小さな値が入っているということはそれだけ速く発散することに相当し、大きな値、特に上限値が入っている場合は実際には
発散せずに収束するとみなせます。つまりマンデルブロ集合の元であるということです。

M = 10000000
maxItr = 255
Mand = Mandelbrot(M, maxItr)

x_lower = -2.5
x_upper = 1.1
y_lower = -1.8
y_upper = 1.8
pixel = 1000
z = [[0 for m in range(pixel+1)] for n in range(pixel+1)]
for i in range(pixel+1):
    re = x_lower + (x_upper - x_lower)/pixel * i 
    for j in range(pixel+1):
        im = y_lower + (y_upper - y_lower)/pixel * j
        z[j][i] = Mand.Calc(re, im)

最後に結果を格納したzをヒートマップとしてmatplotlibで描画してみましょう。
ここではnumpy.meshgridを使ってヒートマップを描いていきます。
以下を参考にしました。

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(pixel+1)
y = np.arange(pixel+1)
Z = np.array(z)
X, Y = np.meshgrid(x, y)

plt.rcParams['figure.figsize'] = 20,20
plt.gca().set_aspect('equal', adjustable='box')
plt.pcolor(X, Y, Z)
plt.colorbar()
plt.hot()
plt.xlabel('Re', fontsize=22)
plt.ylabel('Im', fontsize=22)
plt.xticks(color="None")
plt.yticks(color="None")
plt.tick_params(length=0)
plt.show()

結果

見事にマンデルブロ集合を再現できました。マンデルブロ集合の醍醐味は、拡大した時に現れるフラクタル図形にありますが、それを表現するには手元のPCの計算リソースは少なすぎます。ここでは自らの手で全体描写できたことに満足して、万華鏡のような美しいフラクタル性は動画に譲ります。

おまけ(ジュリア集合)

漸化式はマンデルブロ集合のそれと同じ形のままで、今度は固定の$c \in \mathrm{Z}$を予め定め、$z_0$を変数とした時に漸化式が発散しないような複素平面上の集合はJulia(ジュリア)集合と呼ばれます。これまた発散の速度をヒートマップにすることでグラフを可視化することができ、やはり美しい図形を描くことができます。上の要領がわかれば簡単に描画できるのでここでもやってみます。

class Julia:
    def __init__(self, M, maxItr, c_r, c_i):
        self.M = M
        self.maxItr = maxItr
        self.c_r = c_r
        self.c_i = c_i

    def Calc(self, r, i): # re, im: 2d Array 
        z_r = r
        z_i = i

        for i in range(self.maxItr):
            if z_r**2 + z_i*2 > self.M:
                return i
            z_r_ = z_r**2 - z_i**2 + self.c_r
            z_i_ = 2*z_r*z_i + self.c_i
            z_r, z_i = z_r_, z_i_
        return self.maxItr

Juliaクラスのメンバー変数に、複素定数の$c$をセットできるようにした以外はそれほど大きな変更はありません。

import numpy as np
import matplotlib.pyplot as plt

M = 10000000
maxItr = 255
c_r = 0.06
c_i = -0.70

J = Julia(M, maxItr, c_r, c_i)

x_lower = -2.0
x_upper = 2.0
y_lower = -1.5
y_upper = 1.5
pixel = 1000
z = [[0 for m in range(pixel+1)] for n in range(pixel+1)]
for i in range(pixel+1):
    re = x_lower + (x_upper - x_lower)/pixel * i 
    for j in range(pixel+1):
        im = y_lower + (y_upper - y_lower)/pixel * j
        z[j][i] = J.Calc(re, im)

x = np.arange(pixel+1)
y = np.arange(pixel+1)
Z = np.array(z)
X, Y = np.meshgrid(x, y)

plt.rcParams['figure.figsize'] = 20,20
plt.gca().set_aspect('equal', adjustable='box')
plt.pcolor(X, Y, Z)
plt.hot()
plt.xlabel('Re', fontsize=22)
plt.ylabel('Im', fontsize=22)
plt.xticks(color="None")
plt.yticks(color="None")
plt.tick_params(length=0)
plt.show()

実行してみて、どんな図形が描けるか確認してみてください。特に複素定数の$c$(プログラム中のc_rとc_i)を弄ると、図形の様子が様々に変わると思います。ぜひ美しい図形を探索してみてください。(この場合はポインセチア畑みたいな?のが描けました)