【はじめてのPython】マンデルブロ集合を描く - 出世はあきらめて30代から手に職をつけてみる

こんにちは、ヒトツメです。
今日はPythonの話題なのですが、前回のミニマックス法→アルファ・ベータ法の話はいったん横において、今日は少し違った話題を取り上げようと思います。

何を取り上げるかというと、ブノワ・マンデルブロという数学者が発見した、マンデルブロ集合と呼ばれるものです。
一見ちょっと奇妙な集合ですが、知れば知るほどとりこになる、少し不思議な集合の話です。

複素平面上の集合体
発散しない点の集合
Pythonで実装してみる。
さいごに

複素平面上の集合体

まずはPythonの実装の前に、マンデルブロ集合とは何か？という話ですが、マンデルブロ集合を知るためには、複素平面について知っておかなければなりません。
高校数学で少しかじったという人も多いと思いますが、複素平面とは、小学生の時に習う数直線を平面に拡大したものです。数直線というと、真ん中に「0」、右に行けば行くほどプラス、左に行けば行くほどマイナスになる直線を言います。

19世紀、この数直線には、x方向の軸だけではなく、y方向の軸もあるということがガウスによって提唱されました。そのような考えにしがたい、2乗すると「-1」になる数字「i」をy方向にとることで、数直線を平面に拡大したものを、複素平面と言います。

例えば複素数「2+i」はx軸が2、y軸が1の場所にプロットされることとなります。

発散しない点の集合

さて、そんな複素平面上にある点ですが、ひたすら2乗して元の複素数を足すという作業を繰り返してあげると、無限に近い数字に発散する場合としない場合が存在します。マンデルブロ集合は、この時発散しないものを集めた点の集合を指します。
つまり、次のような式で表される数式において、どのようなCを取れば、nがどのような数字の時に無限に近い数字まで発散してしまうのかをプロットしていったのが、マンデルブロ集合です。

これをヒートマップの形式にして表すと、次のような集合となります。

この図形の面白いところは、この「仏さまが座ったような形の図形」が、色々な場所に出てくるところです。実際、頭の少し上の方（左の方）にとても小さな茶色い点が見えると思いますが、これを大きく拡大すると、ほぼこれと同じ図形が出てきます。
このように自己相似性を持つ図形をフラクタル構造といい、一般的には、野菜のロマネスクの形状がこれに該当することが知られています。

Pythonで実装してみる。

さて、一見するとこのような奇妙な集合も、Pythonでは複素数の計算が簡単にできるので、意外とあっさりと実装することができます。
今回はヒートマップの生成のために、pyplotというライブラリを使いますが、それ以外は、一般によく使うnumpyだけで計算できてしまいます。

import numpy as np
import matplotlib.pyplot as plt

まずは、上記の数式で、任意のCを設定したときに、nがいくつのときに発散するかを計算する関数を導入します。このとき、「発散しない場合にnがいくつの時に計算をやめるか」と、「0からの距離がいくつになったら発散したとみなすか」という二つの問題を事前に決めておく必要があります。
二つ目に関しては、0からの距離が2を超えれば戻ってくることはないということが分かっているので、2を設定します。一つ目の問題については、任意に設定ができるように、関数の中の引数に設定します。

def func_mandelbrot(max, re_num, im_num):
    c = complex(re_num, im_num)
    z = complex(0, 0)
    for i in range(max):
        z = z ** 2 + c
        if abs(z) >= 2:
            return i
    return max

re_numは実数部を、im_numは虚数部を示します。complex関数でCを設定し、Cに対して、引数maxの回数を最大として「2乗してCを足す」という作業を繰り返していきます。この関数を用いることで、何回で発散したかを返却することができます。

あとはこの関数を、事前に指定した範囲における任意のCに対して当てはめていき、それをベースにヒートマップを設定すれば、上記のような図形が出力できます。

def set_mandelbrot(left, right, bottom, top, resolution, max):
    if (right - left) > (top - bottom):
        width = resolution
        height = round(resolution / (right - left) * (top - bottom))
    elif (right - left) < (top - bottom):
        width = round(resolution / (top - bottom) * (right - left))
        height = resolution
    else:
        width = resolution
        height = resolution
    #返却用の配列
    mandelbrot_array = np.zeros((height, width))
    #縦横を解像度に分解
    re_array = np.linspace(left, right, width)
    im_array = np.linspace(top, bottom, height)
    #縦横のポイントごとにfunc_mandelbrotを適用
    for w_num, re_num in enumerate(re_array):
        for h_num, im_num in enumerate(im_array):
            mandelbrot_array[h_num][w_num] = func_mandelbrot(max, re_num, im_num)
    return mandelbrot_array

これで、左端・右端・上部・下部の数字を順に設定すれば、それをresolutionで設定した数字分に分割して、nがどの数字を取る場合に発散するかを示す数式が返却可能です。
ちなみに、一つ工夫をしてあって、縦横のサイズが違う可能性を考慮し、長い方の辺はresolutionの数字で分割しつつ、短い方の辺はその解像度と同じ解像度になるように、分割の数字を減らすように設定しています。

top = 1.2
bottom = -1.2
left = -2.4
right = 0.8
resolution = 800
max = 200
Mandelbrot = set_mandelbrot(left, right, bottom, top, resolution, max)
plt.imshow(Mandelbrot, cmap="jet")
plt.show()

このようにして出力した配列を、pyplotのimshowに取り込めば、集合が出力されます。
ちなみに、”jet”の部分を色々変えると、様々な色合いが設定できます。