Matlab で マンデルブロ集合

息抜きで、Matlabでマンデルブロ集合の図形を描写するプログラムを作成してみました。

マンデルブロ集合とは

マンデルブロ集合は、Wikipediaから以下です。

次の漸化式で定義される複素数列\{z_n\}_{n \in \mathrm{N}}がn\to\inftyの極限で無限大に発散しないという条件を満たす複素数c全体が作る集合がマンデルブロ集合である。
（Wikipediaより引用）

    \[ \left\{ \begin{aligned} z_{n+1} &= {z_n}^2 + c \\ z_0 &= 0 \end{aligned} \right. \]

ざっくりとしたプログラムの解説としては、複素平面座標にて、その複素数の発散の度合いで図形が描ける、というものです。

ソースコード

実際のソースコードは以下です。

複素平面にて描写範囲を指定し、各座標の複素数Cの発散を調べています。
絶対値が2以上になった場合を発散しているとし、その発散回数を保存しています。
図形の描写では、各座標の発散までの計算回数をプロットしています。
（変数は単精度で計算しています）

% ----------------------------------------------------- %
% パラメータ設定
% ----------------------------------------------------- %
% 最大ループ回数
loop = 1e2;
% 各軸の描写点数
n    = 1e3;
% 描写範囲
xmin = -1.5;
xmax = 0.5;
ymin = -1.0;
ymax = 1.0;
% ----------------------------------------------------- %

% ----------------------------------------------------- %
% マンデルブロ集合の計算
% ----------------------------------------------------- %
% 座標設定
x = single(linspace(xmin,xmax,n));
y = single(linspace(ymin,ymax,n));
C = ones(n,1)*x + 1i*y(:)*ones(1,n);
% マンデルブロ図形の計算
Z = zeros(n,'single');
W = zeros(n,'single');
% インデックスの初期値
index = true(size(C));
% ループ
for k = 1:loop
    Z(index) = Z(index).^2 + C(index);
    index = abs(Z) < 2;
    W(~index & W==0) = k;
end
W(W==0) = loop+1;
% 図形のプロット
imagesc(x,y,W);
axis square;
title('Mandelbrot Set');
xlabel('Real');
ylabel('Imag');
% PNG 画像の保存
print('Mandelbrot.png', '-dpng');
% ----------------------------------------------------- %

このソースコードで描いてみた図形が以下。

Mandelbrot1

描写位置を変更すると面白いです。

Mandelbrot2

Mandelbrot3