Processing マンデルブロ集合

こんにちはNUMです。

7月に入ってからとても暑くなってきましたね。僕は最近ZINE制作に夢中になっています。 ZINEとは様々なトピックについて情報共有したり、作家であればアートブックとして販売したりする小雑誌です。 最近行ったON READINGという本屋さんにZINEが沢山置いてあって、色んな作家さんのZINEを見てると、どれも個性的な世界観で本当に見入ってしまうものばかりでした。それに本の品揃え、空間どれも素晴らしいので、ぜひ名古屋にいらした際は行ってみて下さい。

ON READINGの入口

センスの良い本が沢山

今回はジェネラティブアートをやっている人であれば大体知っているであろうマンデルブロ集合を作っていきます。

目次

マンデルブロ集合とは

歴史

プログラムにおける複素数の扱い

概要

実装

まとめ

参考

マンデルブロ集合とは

マンデルブロ集合とは、上記漸化式で表される複素数列 Zn の絶対値 |Zn| がn→∞の極限で無限大に発散しないような複素数Caの集合です。

マンデルブロ集合

数学に親しくない人は難しいと思いますが、要はCに値を当てはめてみてZnが無限大に発散しない集合という理解で良いです。それを可視化すると複雑でありつつ綺麗な模様が浮かび上がるのです。

歴史

マンデルブロ集合の由来になっている数学者ブノア・マンデルブロは「フラクタル」という用語を初めて提唱し、自然界に広く見られる不規則形状を数理的に記述するフラクタル理論を体系化しました。 雲や山岳、樹木の枝分かれなどを解析対象として提示し、視覚的直感を重視した解説で一般読者にも自然と数学の関連性について大きな影響を与えました。 マンデルブロ集合は、20世紀初頭にフランスの数学者ピエール・ファトゥとガストン・ジュリアによって初めて研究された複素力学系*1 に起源を持ちます。 *1 複素数を変数とする関数を繰り返し適用したときの振る舞いを研究する分野 マンデルブロ集合の名前の元となる人物であるブノア・マンデルブロは、1980年に発表された論文で、二次多項式のパラメータ空間を研究しました。実は、マンデルブロ集合の数学的研究はブノア・マンデルブロが最初の研究ではなく、数学者アドリアン・ドゥアディとジョン・H・ハバード（1985年）の研究から始まり、 二人はフラクタル幾何学における影響力のある研究に敬意を表して、この集合にマンデルブロの名称を付けました。

プログラムにおける複素数の扱い

Processing(Java)では複素数を扱うことができないため、実部と虚部に分けた数値演算に落とし込む必要があります。

複素数は実数では表せない虚数を含む数で以下形式で表すことができます。 虚数単位 i は数学で使用する変数と同じ扱いで計算可能です。ただし i 同士をかけると-1になるので注意して下さい。

マンデルブロ集合の数式は全て複素数なので、まずはプログラムに置き換えるためにそれぞれa+bi形式に変換して実部と虚部に分けます。

この実部と虚部をfloat形式でそれぞれ計算していけばマンデルブロ集合は出来上がります。

概要

Processingにおけるマンデルブロ集合は以下方針で実装していきます。

1. キャンバスを格子状に分割

2. 各セル毎に以下処理をループ

   1. 各セルに対応する定数Cの実部、虚部を-2~2の範囲で取得

   2. 手順aで取得した値、前項の値を使用して数列が発散するか判定

3. 発散しない場合は黒、発散する場合はループ回数に応じて白から黒のグラデーションで塗る

   
   

   

   手順a：セルに定数の実部、虚部をマッピング

実装

int N = 200; // ループ回数

void setup() {
  size(850, 850);
  pixelDensity(displayDensity());
  background(255);

  noStroke();
  for (float i = 0; i < width; i++) {
    for (float j = 0; j < height; j++) {
      int result = isDivergence(i, j);
      if (result == 1) {
        fill(0); // // 発散しない（マンデルブロ集合に属する）
      } else {
        fill(map(-result, 0, N, 255, 0)); // ループ回数に応じてグレー値をマッピング
      }
      rect(i, j, 1, 1);
    }
  }
}

/**
 * マンデルブロ集合の発散判定
 * 複素数c = a + ibに対して、Zn+1 = Zn^2 + cをループし
 * どのタイミングで |z| >= 2 となるか判定
 *  
 * @param1　float x座標
 * @param2　float y座標
 * @return　int 1：発散しない 1以外：発散する
 */
int isDivergence(float x, float y) {
  // 初項は0で次項の計算値で更新
  float startRe = 0; // 初項の実部
  float startIm = 0; // 初項の虚部

  // 定数の実部cと虚部dを-2~2の範囲に座標にマッピング
  float c = map(x, 0, width, -2, 2);
  float d = map(y, 0, height, -2, 2);

  for (int i = 0; i < N; i++) {
    // マンデルプロ集合の計算
    float re = calcRealPart(startRe, startIm, c); // 実部
    float im = calcImaginaryPart(startRe, startIm, d); // 虚部

    // 複素数の大きさを計算
    float mag = sqrt((re * re) + (im * im));

    // 発散の判定 zが2以上であれば発散
    if (mag >= 2) {
      return -i; // fillの色を変化させるために、何回目のループで発散したかを返す
    }
    // 次項の計算のために更新
    startRe = re;
    startIm = im;
  }  
  
  // 最大ループまで発散しなかった場合、マンデルブロ集合に属する
  return 1;
}

/**
 * 実部の計算
 * @param1　float 計算項の実部
 * @param2　float 計算項の虚部
 * @param3　float 定数項の実部
 * return   float 実部の計算値
 */
float calcRealPart(float re, float im, float c) {
  return (re * re) - (im * im) + c;
}

/**
 * 虚部の計算
 * @param1　float 計算項の実部
 * @param2　float 計算項の虚部
 * @param3　float 定数項の虚部部
 * return   float 虚部の計算値
 */
float calcImaginaryPart(float re, float im, float d) {
  return (2 * re * im) + d;
}

それではプログラムを理解する上で重要なポイントに絞って解説していきます。

for (float i = 0; i < width; i++) {
  for (float j = 0; j < height; j++) {
    int result = isDivergence(i, j);
    if (result == 1) {
      fill(0); // // 発散しない（マンデルブロ集合に属する）
    } else {
      fill(map(-result, 0, N, 255, 0)); // ループ回数に応じてグレー値をマッピング
    }
    rect(i, j, 1, 1);
  }
}

キャンバスを格子状に分割し、各セルでマンデルブロ集合の数式が発散するか判定をします。 発散する場合はセルを黒色、発散しない場合はループ回数に応じて黒から白へグラデーションさせていきます。

isDivergence関数の戻り値がループ回数をマイナスの値で返すので、マイナス符号をresultに付与して符号を反転させてグラデーションの値に使用します。

/**
 * マンデルブロ集合の発散判定
 * 複素数c = a + ibに対して、Zn+1 = Zn^2 + cをループし
 * どのタイミングで |z| >= 2 となるか判定
 *  
 * @param1　float x座標
 * @param2　float y座標
 * @return　int 1：発散しない 1以外：発散する
 */
int isDivergence(float x, float y) {
  // 初項は0で次項の計算値で更新
  float startRe = 0; // 初項の実部
  float startIm = 0; // 初項の虚部

  // 定数の実部cと虚部dを-2~2の範囲に座標にマッピング
  float c = map(x, 0, width, -2, 2);
  float d = map(y, 0, height, -2, 2);

  for (int i = 0; i < N; i++) {
    // マンデルプロ集合の計算
    float re = calcRealPart(startRe, startIm, c); // 実部
    float im = calcImaginaryPart(startRe, startIm, d); // 虚部

    // 複素数の大きさを計算
    float mag = sqrt((re * re) + (im * im));

    // 発散の判定 zが2以上であれば発散
    if (mag >= 2) {
      return -i; // fillの色を変化させるために、何回目のループで発散したかを返す
    }
    // 次項の計算のために更新
    startRe = re;
    startIm = im;
  }  
  
  // 最大ループまで発散しなかった場合、マンデルブロ集合に属する
  return 1;
}

isDivergence関数はマンデルブロ集合が発散するかを判定する中核となる関数です。

セルの座標から定数Cの実部、虚部を取得する方法は以下です。

// 定数の実部cと虚部dを-2~2の範囲に座標にマッピング
float c = map(x, 0, width, -2, 2);
float d = map(y, 0, height, -2, 2);

数列の計算は以下関数で実部、虚部をそれぞれ実施しています。

/**
 * 実部の計算
 * @param1　float 計算項の実部
 * @param2　float 計算項の虚部
 * @param3　float 定数項の実部
 * return   float 実部の計算値
 */
float calcRealPart(float re, float im, float c) {
  return (re * re) - (im * im) + c;
}

/**
 * 虚部の計算
 * @param1　float 計算項の実部
 * @param2　float 計算項の虚部
 * @param3　float 定数項の虚部部
 * return   float 虚部の計算値
 */
float calcImaginaryPart(float re, float im, float d) {
  return (2 * re * im) + d;
}

複素数の大きさは三平方の定理で計算し、2以上であれば発散したとみなし、その時点のループ回数をマイナス符号を付与して返却します。

マイナス符号を付与する理由は、発散しない場合に1を返すので、たまたま2回目(ループは0回からなので)のループで発散した場合に戻り値が同じになるため、誤判定を防ぐためです。

// 複素数の大きさを計算
float mag = sqrt((re * re) + (im * im));

// 発散の判定 zが2以上であれば発散
if (mag >= 2) {
  return -i; // fillの色を変化させるために、何回目のループで発散したかを返す
}

まとめ

今回はマンデルブロ集合の解説記事でした。実は、マンデルブロ集合は定数の範囲を調整して描画する部分を拡大すると様々な模様が現れることが分かっています。

cとdを以下範囲にすると以下のような模様が出来上がります。（色は違いますが）

c = -0.36 ~ -0.37;

d = -0.64 ~ -0.65;

パラメータを調整するだけで全然違った見え方になるのでずっと遊んでいられます。 ぜひ皆さんも自分だけのパラメータを見つけて下さい！

最後まで読んでいただきありがとうございました！

参考

• https://handwiki.org/wiki/Mandelbrot_set#History

• https://azisava.sakura.ne.jp/mandelbrot/definition.html