マンデルブロ集合 - マンデルブログ

やったったぜ

昨日、やると言ってやらなかったマンデルブロ集合のプログラム
今日こそ作りますよ。

マンデルブロ集合とは

これがマンデルブロ集合です。
一応簡単に説明をしておきますと、
まず $z_0=0$ として漸化式 $z_n=z_{n+1} +c$ を考えます。 $c$ は複素数です。
$n\to\infty$ の極限を取った時、 $z_n$ が無限大に発散しない $c$ の集合がマンデルブロ集合です。
この複素数 $c$ を複素平面上の点 $c=a+ib$ だと思ってプロットするとフラクタル集合になっています。

あんまり説明できてない気がするので、詳しくはWikipediaを見てください。

いざプログラミング

いつもはatomで書いてるんですが、
お家のパソコンで書いたのでJupyter Notebookを使いました。
以下のブログを参考にさせていただきました。
Pythonでカオス・フラクタルを見よう！

import numpy as np
import matplotlib.pyplot as plt

def mandelbrot(array):
    a = array[0]
    b = array[1]
    x = 0
    y = 0
    N = 50
    for i in range(N):
        x,y = x**2 - y**2 + a , 2*x*y + b
        try:
            if x**2 + y**2 > 4:
                return 1 - (i/N)
        except:
            break
        
    return 0

$z_{n}=x_{n}+iy_{n}$ , $c=a+ib$ としたとき、
$z_{n+1}^2=(x_{n}+iy_{n})(x_{n}-iy_{n})=(x_{n}^2-y_{n}^2)+i(2x_{n}y_{n})$
よって
$z_{n+1}^2+c=(x_{n}^2-y_{n}^2+a)+i(2x_{n}y_{n}+b)$
となり、
$x_{n+1}=x_{n}^2-y_{n}^2+a$
$y_{n+1}=2x_{n}y_{n}+b$
となります。数式を中央揃えにするのってどうやるんだろ。

この関数は点の座標を代入すると、発散具合を返してくれる関数です。
速く発散するほど1に近い値になり、発散しない点では0を返します。
無限回繰り返すわけにはいかないので、参考にさせていただいたページ通り
適当な回数繰り返して $x^2 + y^2 >4$ であれば発散としています。

x = [np.linspace(-2, 1, 500)]
y = [np.linspace(-1.5,1.5,500)]
X, Y = np.meshgrid(x, y)

集合の全ての座標を書き出します。
numpy.meshgrid()は代入した2つの配列の全ての組み合わせを返してくれる便利な関数です。

mesh_arr = np.c_[X.ravel(),Y.ravel()]
color = [mandelbrot(mesh_arr[i]) for i in range(len(mesh_arr))]
color= np.reshape(color,[len(X),len(Y)])

mesh_arrで[x,y]を縦に並べています。
colorはその座標の色(発散具合）をため込むlistです。
numpy.reshape()で二次の形に直しています。

plt.figure(figsize=(8,6))
plt.pcolor(X, Y, color)
plt.savefig("mandelbrot_set.png",bbox_inches="tight", pad_inches=0.0)

描画して保存です。
こうしてできたマンデルブロ集合がこちらになります。

f:id:woodhero0908:20180606233106p:plain — 色が濃い部分がマンデルブロ集合。

おー、短時間で作った割には結構綺麗にできましたね。
もっと簡単で効率がいい方法もあると思いますが、今回は満足。