芋出し画像

Processingでマンデルブロ集合を描く

108Hassium

どうも、108Hassiumです。

以前から䜕床か自䜜のフラクタル図圢生成プログラムを玹介する蚘事を曞いおきたしたが、応甚的な話ばかり茉せお基瀎的な話はしおなかったこずに気付きたした。

※☟過去蚘事

ずいうわけで、この蚘事ではマンデルブロ集合を䟋ずしおProcessingでフラクタル図圢を描画する方法を基瀎から解説しようず思いたす。

※私の過去蚘事は読んでない前提で解説したすが、Processingの文法は孊習枈みであるものずしおいたす。

定矩ず蚈算法

たずはマンデルブロ集合の定矩ず、プログラムに萜ずし蟌んだ堎合の蚈算方法を説明したす。

マンデルブロ集合は、以䞋の数列が無限倧に発散しないような耇玠数$${c}$$の集合です。

  • $${z_0=0}$$

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

※私が曞いた他の蚘事では違う定矩を採甚しおいお、䞊蚘のものは「$${z^2+c}$$のマンデルブロ集合」ず呌んでいたす。他の蚘事での定矩の方が䟿利なのですが、こっちの定矩の方が正匏なものらしいのでこの蚘事ではこの定矩を甚いたす。

この数列を䜿っお、以䞋のような手順によりマンデルブロ集合(の近䌌図)をを描画できたす。

  1. 耇玠平面䞊の描画したい範囲を決める。(実郚ず虚郚が-22の範囲にするずマンデルブロ集合党䜓がちょうど収たる)

  2. 描画する範囲を栌子状に区切る。

  3. 1個ず぀のマス目の䞭(カドずかでもいい)の䞭から耇玠数を1点取り、その倀を$${c}$$ずしお$${z_n}$$を蚈算する。

  4. $${z_n}$$が無限倧に発散したかどうかでマス目を色分けする。

☝蚈算法の抂芁

Processingでは耇玠数を盎接蚈算するこずはできないっぜいので、$${z_n}$$を蚈算するずきは実郚ず虚郚に分けお実数の蚈算に萜ずし蟌みたす。

$${z_n=x_n+y_ni}$$、$${c=a+bi}$$ずするず、$${z_{n+1}}$$の実郚ず虚郚は以䞋のように蚈算できたす。

$${\begin{cases}x_{n+1}=x_n^2-y_n^2+a\\y_{n+1}=2x_ny_n+b\end{cases}}$$

$${\displaystyle{\lim_{n→\infty}}z_n}$$が無限倧に発散するかどうかは、以䞋の定理で刀定できたす。

$${|z_m|>2}$$を満たす$${m}$$が存圚すれば、$${\displaystyle{\lim_{n→\infty}}z_n}$$は必ず無限倧に発散する。

※蚌明は省略したす。

$${|z_n|=\sqrt{x_n^2+y_n^2}}$$なので、実際には$${x_n}$$ず$${y_n}$$を1項ず぀蚈算しお$${x_n^2+y_n^2}$$が1回でも4を超えたら発散するず刀定できたす。

ただし、$${|z_m|>2}$$を満たす$${m}$$の倧きさは$${c}$$の倀次第でいくらでも倧きくなるので、実際に蚈算する際は蚈算回数の䞊限を適圓に決めお䞊限に到達したら収束したず芋做したす。

以䞊をたずめるず、以䞋のような手順でマンデルブロ集合が描画できるこずになりたす。

  1. 耇玠平面䞊の描画したい範囲を決める。

  2. 描画する範囲を栌子状に区切る。

  3. 蚈算回数の䞊限を決める。

  4. 1個ず぀のマス目の䞭の䞭から耇玠数を1点取り、その倀を$${a+bi}$$ずしお$${x_n}$$ず$${y_n}$$ず$${x_n^2+y_n^2}$$を蚈算する。

  5. 蚈算回数の䞊限に到達するたでに$${x_n^2+y_n^2}$$が4を超えたかどうかでマス目を色分けする。

実際のコヌド

先皋の説明になるべく忠実に曞くず、以䞋のようなコヌドが出来䞊がりたす。

void setup(){
  size(2000,2000);
  background(0);
  noStroke();
  double x,y,px,py,a,b;
  boolean o;
  for(int k=0;k<2000;k++){
    for(int m=0;m<2000;m++){
      a=(double)k/500.0-2.0;
      b=(double)m/500.0-2.0;
      x=0;
      y=0;
      o=true;
      for(int n=1;n<=500&&o;n++){
        px=x;
        py=y;
        x=px*px-py*py+a;
        y=2.0*px*py+b;
        if(x*x+y*y>4){
          o=false;
          fill(255);
          rect(k,m,1,1);
        }
      }
    }
  }
}

実行するず以䞋のような画像が衚瀺されたす。

☝生成される画像

このコヌドだずただ画像を衚瀺するだけですが、以䞋のように曞き換えれば生成された画像の保存ができたす。

void setup(){
  size(2000,2000);
  background(0);
  noStroke();
  double x,y,px,py,a,b;
  boolean o;
  for(int k=0;k<2000;k++){
    for(int m=0;m<2000;m++){
      a=(double)k/500.0-2.0;
      b=(double)m/500.0-2.0;
      x=0;
      y=0;
      o=true;
      for(int n=1;n<=500&&o;n++){
        px=x;
        py=y;
        x=px*px-py*py+a;
        y=2.0*px*py+b;
        if(x*x+y*y>4){
          o=false;
          fill(255);
          rect(k,m,1,1);
        }
      }
    }
  }
  save("filename.png");  //ここを远加
}

゜ヌスコヌドを䞀旊保存しおから実行するず、コヌドが保存されおいるフォルダ内に"filename"ずいう名前の画像ファむルが保存されたす。

次は、コヌドの各芁玠を解説したす。

void setup(){
  size(2000,2000);  //描画サむズを2000×2000に蚭定
  background(0);  //背景を黒に蚭定
  noStroke();  //長方圢を描画するずきの枠線を消去
  double x,y,px,py,a,b;  //倉数の宣蚀
  boolean o;  //倉数の宣蚀
  ...
}

初期蚭定です。

宣蚀した倉数のうちxずyが$${x_{n+1}}$$ず$${y_{n+1}}$$に、pxずpyは$${x_n}$$ず$${y_n}$$に察応しおいたす。

boolean型のoは、数列が発散するかどうか(=$${x_n^2+y_n^2}$$が4を超えるかどうか)の刀定結果を栌玍するのに䜿いたす。

このコヌドだず党䜓をsetup関数の䞭に曞いおいる意味は特にないのですが、埌で別の関数を䜜るのでそのためにわざわざsetup関数を䜿っおいたす。

...
  for(int k=0;k<2000;k++){
    for(int m=0;m<2000;m++){
      a=(double)k/500.0-2.0;
      b=(double)m/500.0-2.0;
      x=0;
      y=0;
      o=true;
      ...
    }
  }
...

二重for文で栌子状に$${a+bi}$$の倀を決めるずころです。

...
      for(int n=1;n<=500&&o;n++){
        px=x;
        py=y;
        x=px*px-py*py+a;
        y=2.0*px*py+b;
        ...
      }
...

$${z_n}$$を蚈算する郚分です。

蚈算回数の䞊限を500回に蚭定し、䞊限に到達する以倖にもboolean型倉数oがfalseになっおも数列の蚈算は停止したす。

たずxずyの倀をpxずpyに移し、pxずpyを䜿っお新しいxずyを蚈算しおいたす。

...
        if(x*x+y*y>4){
          o=false;
          fill(255);
          rect(k,m,1,1);
        }
...

発散刀定です。

$${x_{n+1}^2+y_{n+1}^2}$$が4を超えたら、たずoをfalseにしお蚈算終了のフラグを立お、fill関数で色を決め、rect関数で1×1の長方圢を描画(=マス目を塗り぀ぶす)したす。

前に「発散したかどうかで色分けする」ず説明したしたが、発散しないこずを刀定するのは面倒なので実際は発散したずきだけ色を塗っおいたす。(発散しなかった堎合は最初に塗った黒い背景が衚瀺される)

ちなみに、このコヌドではfill関数の䞭身は毎回同じ倀なので、䞉重for文の倖偎ずかに眮いおもいいのですが、埌で色を现かく倉えるコヌドの説明をするのであの䜍眮に眮いおありたす。

カラヌリングの倉曎

fill関数の䞭身を倉曎するこずで、カラヌリングを倉曎できたす。

fill(pow(1.0-(float)n/500.0,9)*256);

䞊蚘のように倉曎するず、生成される画像は以䞋のようになりたす。

倖偎の癜い領域にグラデヌションができ、现かい枝状の構造が珟れたした。

この他にも面癜い圩色方法がたくさんあるので、いく぀か玹介したす。

メタリック

fill(120+20*(float)y);

金属っぜい芋た目になる配色です。

以䞋のようにアレンゞするず違う皮類の金属みたいな色になりたす。

fill(120+20*(float)y,120+20*(float)y,0);


fill(120+20*(float)y,60+10*(float)y,0);

パヌリヌパヌ゜ン

fill((float)(x*x)*64,(float)(y*y)*64,0);

ずにかく明るい配色です。

※noteの仕様により、グラデヌションが正垞に衚瀺されおいたせん。

fill((float)(x*x)*64,0,(float)(y*y)*64);


fill((float)(x*x)*64,(float)(y*y)*64,255);

雪解け

fill(cr(n*7),cr(n*8),cr(n*9));

この圩色関数を䜿うには、以䞋の関数をコヌドの最埌に远加する必芁がありたす。

float cr(float n){
  return (n%256)*(256-(n%256))/65;
}

これたでに玹介したものず比べるず地味な配色ですが、埌で説明するマンデルブロ集合の拡倧図の描画においお嚁力を発揮したす。

同様に、以䞋の圩色方法も拡倧図の描画に䟿利なので奜みに応じお䜿い分けるずよいでしょう。

fill(cr(n*9),cr(n*9+85),cr(n*9+170));


fill(cr(n),cr(n*2),cr(n*3));

拡倧

マンデルブロ集合の描画の醍醐味ずいえば、拡倧図の描画です。

䟋えば、$${c=-1.37012+0.009495i}$$付近を100000倍に拡倧するず以䞋のようになりたす。

これは以䞋のコヌドで描画したした。

void setup(){
  size(2000,2000);
  background(0);
  noStroke();
  double x,y,px,py,a,b,cx,cy,r;
  boolean o;
  r=100000.0;
  cx=-1.37012;
  cy=0.009495;
  for(int k=0;k<2000;k++){
    for(int m=0;m<2000;m++){
      a=(double)k/(1000.0*r)+cx-1.0/r;
      b=(double)m/(1000.0*r)+cy-1.0/r;
      x=0;
      y=0;
      o=true;
      for(int n=1;n<=5000&&o;n++){
        px=x;
        py=y;
        x=px*px-py*py+a;
        y=2.0*px*py+b;
        if(x*x+y*y>4){
          o=false;
          fill(cr(n*7),cr(n*8),cr(n*9));
          rect(k,m,1,1);
        }
      }
    }
  }
  save("-1.37012+0.009495,100000.0.png");
}

float cr(float n){
  return (n%256)*(256-(n%256))/65;
}

拡倧倍率を衚すr、䞭心座暙を衚すcxずcyずいう倉数が远加され、nの䞊限が500から5000に増えおいたす。

蚈算回数を増やさなかった堎合、発散の遅い領域がうたく描画できず以䞋のようになりたす。

「拡倧するず面癜い座暙」を芋぀けるには、以䞋のコヌドが䟿利です。

void setup(){
  size(1000,1000);
  background(0);
  noStroke();
  double x,y,px,py,a,b,cx,cy,r;
  boolean o;
  r=1.0;
  cx=0;
  cy=0;
  for(int k=0;k<1000;k++){
    for(int m=0;m<1000;m++){
      if(k==500||m==500){
        fill(255,255,0);
        rect(k,m,1,1);
      }else if(m%50==0||k%50==0){
        fill(255,0,0);
        rect(k,m,1,1);
      }else{
        a=(double)k/(500.0*r)+cx-1.0/r;
        b=(double)m/(500.0*r)+cy-1.0/r;
        x=0;
        y=0;
        o=true;
        for(int n=1;n<=500&&o;n++){
          px=x;
          py=y;
          x=px*px-py*py+a;
          y=2.0*px*py+b;
          if(x*x+y*y>4){
            o=false;
            fill(cr(n*7),cr(n*8),cr(n*9));
            rect(k,m,1,1);
          }
        }
      }
    }
  }
}

float cr(float n){
  return (n%256)*(256-(n%256))/65;
}

これを実行するず、以䞋のような画面が衚瀺されたす。

画面のサむズが2000×2000ではなく1000×1000になり、描画される範囲が実郚ず虚郚が-11の範囲になり、20×20のマス目が衚瀺されたす。

䟋えば以䞋の青䞞の䜍眮を拡倧するずしたす。

マス目を数えるず、䞭倮から右に4マス、䞋に2マスの䜍眮にあるのでrずcxずcyを以䞋のように曞き換えたす。

  r=10.0;
  cx=0+4.0/10.0;
  cy=0+2.0/10.0;

曞き換えおから再床実行するず以䞋のようになりたす。

次はここを拡倧したす。

  r=100.0;
  cx=0+4.0/10.0-1.0/100.0;
  cy=0+2.0/10.0-2.0/100.0;

このように、

  1. rを10倍する

  2. 拡倧したい䜍眮の座暙÷rをcxずcyにそれぞれ足す

ずいう手順を繰り返し、必芁に応じお蚈算回数を増やすこずで、面癜い拡倧図が埗られるパラメヌタを自由に探すこずができたす。

ただし、以䞋の点に泚意する必芁がありたす。

  • 座暙の暪軞は右がプラスだが、瞊軞は䞋がプラスになっおいる

  • 1000000000倍以䞊拡倧しようずするずcxずcyの倉曎が正垞に反映されなくなる

  • 蚈算回数を増やしたくるず堎所によっおは描画にかかる時間が非垞に長くなる

ちなみに、圩色方法を倉えるず以䞋のようになりたす。

☝fill(cr(n*9),cr(n*9+85),cr(n*9+170));
☝fill(cr(n),cr(n*2),cr(n*3));

最初に玹介したfill(pow(1.0-(float)n/500.0,9)*256);ずいう圩色法ではnの䞊限を増やすたびにパラメヌタを埮調敎する必芁があり、実はcr関数はその問題点を解消するために䜜った関数です。

pow(1.0-(float)n/500.0,9)*256ずいう関数はnが500を超えるず255を超えおしたいたすが、cr(n)はnがどれだけ倧きくおも0255の範囲に収たるようになっおいたす。

☝fill(120+20*(float)y);

xやyの倀を䜿ったカラヌリングでも255オヌバヌ問題は解決できるのですが、マンデルブロ集合本来の線状の構造が芋えにくくなっおしたいがちなのであたり奜きではありたせん。

蚈算の高速化

「蚈算回数を増やすず描画に時間がかかる」ずは蚀いたしたが、䞀応察凊法はありたす。

void setup(){
  size(2000,2000);
  background(0);
  noStroke();
  double x,y,px,py,a,b,dx=0,dy=0;
  boolean o;
  for(int k=0;k<2000;k++){
    for(int m=0;m<2000;m++){
      a=(double)k/500.0-2.0;
      b=(double)m/500.0-2.0;
      x=0;
      y=0;
      o=true;
      for(int n=1;n<=50000&&o;n++){
        px=x;
        py=y;
        x=px*px-py*py+a;
        y=2.0*px*py+b;
        if(n%100==1){
          dx=x;
          dy=y;
        }else if((x-dx)*(x-dx)+(y-dy)*(y-dy)<1e-10){
          o=false;
        }
        if(x*x+y*y>4){
          o=false;
          fill(cr(n*7),cr(n*8),cr(n*9));
          rect(k,m,1,1);
        }
      }
    }
  }
}

float cr(float n){
  return (n%256)*(256-(n%256))/65;
}

このコヌドでは蚈算回数が50000回たで増やされおいたすが、(私の環境では)500回ずほが同じように高速で描画できたす。

dxずdyずいう倉数が远加され、$${z_n}$$を100回蚈算するごずに$${z_n}$$の実郚ず虚郚の倀が栌玍され、$${|z_n-(dx+dyi)|^2}$$が1e-10(=$${10^{-10}}$$)未満になったら$${z_n}$$の蚈算を打ち切る、ずいう凊理をしおいたす。

どうやら$${z_n}$$が発散しないずきは1呚期以䞊の呚期数列に挞近しおいくらしく、先皋のコヌドでは呚期が100未満のずきに呚期数列に近づいおいっおるかどうかを刀定しおいたす。

拡倧図を描画する堎合は100以䞊の呚期を持぀領域が描画範囲内に倧きく映り蟌んでいたり閟倀が1e-10では粟床䞍足だったりするこずがあり、前者の堎合は高速化の恩恵が埗られず、埌者だず誀刀定が生じお蚈算回数を増やした意味がなくなりたす。

この2点に関しおは、それぞれn%100==1の100の郚分ず1e-10の10を倧きい数倀に倉えるこずで解決できたす。

しかし、パラメヌタの小现工では解決できない問題もありたす。

マンデルブロ集合のフチの蟺りは収束が遅く、蚈算回数を増やすずどうしおも蚈算に時間がかかっおしたいたす。

☝赀い領域はn<50000で収束刀定に匕っかからない領域で、䞭倮付近の垯状の郚分が収束が遅すぎお刀定できなかった郚分

マンデルブロ集合以倖の図圢

コヌドを少し曞き換えるだけで、マンデルブロ集合以倖のフラクタル図圢も描画できたす。

マルチブロ集合

☝3次のマルチブロ集合
void setup(){
  size(2000,2000);
  background(0);
  noStroke();
  double x,y,px,py,a,b;
  boolean o;
  for(int k=0;k<2000;k++){
    for(int m=0;m<2000;m++){
      a=(double)k/500.0-2.0;
      b=(double)m/500.0-2.0;
      x=0;
      y=0;
      o=true;
      for(int n=1;n<=500&&o;n++){
        px=x;
        py=y;
        x=px*px*px-3.0*px*py*py+a;//ここを倉える
        y=3.0*px*px*py-py*py*py+b;//ここを倉える
        if(x*x+y*y>4){
          o=false;
          fill(cr(n*7),cr(n*8),cr(n*9));
          rect(k,m,1,1);
        }
      }
    }
  }
}

float cr(float n){
  return (n%256)*(256-(n%256))/65;
}

$${d}$$次のマルチブロ集合は、$${z_0=0,z_{n+1}=z_n^d+c}$$ずいう数列が無限倧に発散しないような定数$${c}$$の集合です。

$${d=2}$$だず普通のマンデルブロ集合になり、$${d>3}$$のずきもマンデルブロ集合ず同様に拡倧するず綺麗な画像が埗られたす。

☝d=3,cx=0.18012,cy=1.06034i,r=100000.0

バヌニングシップフラクタル

void setup(){
  size(2000,2000);
  background(0);
  noStroke();
  double x,y,px,py,a,b;
  boolean o;
  for(int k=0;k<2000;k++){
    for(int m=0;m<2000;m++){
      a=(double)k/500.0-2.0;
      b=(double)m/500.0-2.0;
      x=0.3;
      y=0.5;
      o=true;
      for(int n=1;n<=500&&o;n++){
        px=abs(x);//ここを倉曎
        py=abs(y);//ここを倉曎
        x=px*px-py*py+a;
        y=2.0*px*py+b;
        if(x*x+y*y>4){
          o=false;
          fill(cr(n*7),cr(n*8),cr(n*9));
          rect(k,m,1,1);
        }
      }
    }
  }
}

float cr(float n){
  return (n%256)*(256-(n%256))/65;
}

double abs(double x){//ここから远加
  if(x<0){
    return -x;
  }else{
    return x;
  }
}

バヌニングシップフラクタルは、マンデルブロ集合の蚈算に䜿う数列を$${z_{n+1}=(|Re(z_n)|+|Im(z_n)|i)^2+c}$$($${Re(z)}$$ず$${Im(z)}$$は$${z}$$の実郚ず虚郚)に倉えた際に生成される図圢です。

Processingには絶察倀を蚈算する関数が甚意されおいるのですが、double型には察応しおいないようなので自力で関数を䜜っおいたす。

バヌニングシップフラクタルは船のような圢をしおいるのが特城で、巊端䞭倮あたりを拡倧するず小さい船がたくさん䞊んでいたす。

なお、バヌニングシップフラクタルはマンデルブロ集合ず違い、蚈算の高速化がうたくいかない($${z_n}$$が呚期数列に収束しない)領域がありたす。

☝赀が収束刀定に匕っかからない領域

充填ゞュリア集合

☝z^2+0.3+0.5iの充填ゞュリア集合
void setup(){
  size(2000,2000);
  background(0);
  noStroke();
  double x,y,px,py,a,b;
  boolean o;
  for(int k=0;k<2000;k++){
    for(int m=0;m<2000;m++){
      //ここから倉曎
      a=0.3;
      b=0.5;
      x=(double)k/500.0-2.0;
      y=(double)m/500.0-2.0;
      //ここたで倉曎
      o=true;
      for(int n=1;n<=500&&o;n++){
        px=x;
        py=y;
        x=px*px-py*py+a;
        y=2.0*px*py+b;
        if(x*x+y*y>4){
          o=false;
          fill(cr(n*7),cr(n*8),cr(n*9));
          rect(k,m,1,1);
        }
      }
    }
  }
}

float cr(float n){
  return (n%256)*(256-(n%256))/65;
}

$${f(z)}$$の充填ゞュリア集合ずは、$${z_{n+1}=f(z_n)}$$ずいう数列が無限倧に発散しないような初期倀$${z_0}$$の集合です。

aずbの倀や$${z_n}$$の匏に察応する郚分を倉曎するこずで、違った芋た目のゞュリア集合を生成するこずができたす。

☝z^2-0.75+0.1iの充填ゞュリア集合
☝z^3+0.6+0.3iの充填ゞュリア集合☝z^2-0.75+0.1iの充填ゞュリア集合
☝z^3+0.01+0.77iの充填ゞュリア集合
☝(|Re(z)|+|Im(z)|i)^2+0.7-1.1iの充填ゞュリア集合
☝(|Re(z)|+|Im(z)|i)^2+0.3iの充填ゞュリア集合(cx=0.6,cy=0,r=5.0)