DXライブラリでマンデルブロ集合をかく

11213 ワード
DXライブラリ 数学 C++ C++ テキストリンク

こんにちは．はっちーです．Cライクな僕にとってDXライブラリは非常に貴重なライブラリです．簡単なわりに有能ですからね！
Unityの勉強もそろそろしたいと思ってるんですけど(笑)
今日はそのDXライブラリを使ってマンデルブロ集合を書いてみました．
マンデルブロ集合，プログラマはやたら好きかもしれません(笑)

マンデルブロ集合の定義

以下はマンデルブロ集合の定義です．．

漸化式

$$z_{n+1}=z_n^2+c $$
$$z_1=0$$

で定義される複素数列が$n \to \infty $で発散しないとき複素数$c$がマンデルブロ集合に属するといい，そのようなc全体のことをマンデルブロ集合という．

実装

マンデルブロ集合の定義をプログラムしやすいように変形してみましょう．
$z_n=x_n+iy_n$として，実部と虚部に分ける操作を行います．
すると，漸化式は

\begin{align}
x_n &= x_n^2-y_n^2+Re(c) \\
y_n &= 2x_ny_n +Im(c)\\
x_0 &= y_0 =0
\end{align}

のように，実数のみの演算で行うことができます．
次に，数列の収束条件です．マンデルブロ集合が収束する範囲は$|c|<2$の範囲であることが分かっています．
もし，$|z_n|>2$ならば，$|z_{n}|^2>4$となります．$c>-2$はですから定義より$|z_{n+1}|>2$を導けます．
これにより，この数列は一度$|z_n|>2$となればそれ以降も増大していく．すなわち発散するということです．
$|z_n| = |x_n+iy_n|=\sqrt{x_n^2+y_n^2}$となります．sqrtを計算する時間を省くために両辺二乗します．
まとめると得られた発散の条件は

x_n^2+y_n^2 > 4

となります．この$n$について色分けするとマンデルブロ集合になります．

コンピュータにおいては$n$がある回数以上になれば収束と判定するという妥協をしまければいけません．
数列の一般項が見つかればなんとでもなるのですがこの漸化式において一般項は見つかっていません．（僕も見つけようとして一晩数式をこねくり回していたことがあります．）

以下にある$c$について，マンデルブロ集合に属しているか属していないか判定する関数を示します．

mandelbrot_function.cpp
#define Judge 200

int Mandelbrot(double x,　double y) {
    double x_0 = 0.0;
    double y_0 = 0.0;
    double x_n, y_n;
    int count = 1;
    while (count < JUDGE) {
        x_n = x_0*x_0 - y_0*y_0 + x;
        y_n = 2.0*x_0*y_0 + y;

        if (x_n*x_n + y_n * y_n > 4.0) return count;    //必ず発散する

        x_0 = x_n;　y_0 = y_n;
        count++;
    }
    return 0;
}

大文字で書いてるのはすべて定数です．
色の付け方は適当です．いい色の付け方があったら教えてください．

mandelbrot.cpp
int WINAPI WinMain(HINSTANCE, HINSTANCE, LPSTR, int) {
    ChangeWindowMode(TRUE), SetGraphMode(X_Gr, Y_Gr, 32);
    DxLib_Init(), SetDrawScreen(DX_SCREEN_BACK);

    double x_min = -2.0;
    double x_max = 2.0;
    double y_min = -2.0;
    double y_max = 2.0;

        for (int i = 0; i < X_Gr; i++) {
            for (int j = 0; j < Y_Gr; j++) {
                long double x, y;
                int color;

                ij2xy(&x, &y, i, j, x_min, x_max, y_min, y_max, X_Gr, Y_Gr);

                color = Mandelbrot(x, y);
                if (color == 0) DrawPixel(i, j, GetColor(0, 0, 0));
                else           DrawPixel(i, j, GetColor(color % 256 + 1, (color % 16) * 16 + 1, (color % 32) * 8 + 1));
                if (!Myloop()) { DxLib_End(); return 1; }
            }

            //DrawString(0, 0 , temp_x, GetColor(0, 0, 0));
            //DrawString(0, 20, temp_y, GetColor(0, 0, 0));
            //if(i%20==0) ScreenFlip();
        }
        ScreenFlip();
    SaveDrawScreenToPNG(0, 0, X_Gr, Y_Gr, "mandelbrot.png", 20);
    DxLib_End();
    return 0;
}

とりあえず，このコードでできたマンデルブロ集合はこちら

綺麗ですね！少し解説するとすれば，黒のところが収束してしまうところ色つきは発散しないところです．

長くなるのでコードは書きませんが永遠とズームしていくプログラムも書きました．gif画像です．
$z= -0.75147732885+0.031315348 i$に近づいていくところです．数値は適当です
浮動小数点数の限界を達したので途中で打ち切ってます．

何か質問やアドバイスなどあればよろしくお願いします．