塵芥回顧録

なるべく更新していきたいが、ネタがない。

p5jsでマンデルブロ集合を描こう

p5jsでマンデルブロ集合を描いてみよう

マンデルブロ集合


目次


マンデルブロ集合とは

マンデルブロ集合 - Wikipedia
マンデルブロ集合とはリンク先の定義にあるように、
以下の数列において発散しない複素数cをプロットした図です。
{\displaystyle 
\begin{eqnarray}
  \left\{
    \begin{array}{l}
      z_{n+1} = z_n^2 + c\\
      z_0 = 0
    \end{array}
  \right.
\end{eqnarray}
}
図では発散しない複素数cを黒、それ以外を発散としています。
今回実数をx軸、虚数をy軸としてプロットしていくため、式を変換していく必要がある。
複素数cc=a+ib
z_nの実部をx_n、虚部をy_nとすると、
{\begin{align}
z_1 &=a+ib
\end{align}
}
{\begin{align}
z_2 &=(a+ib)^2+a+ib\\
&=a^2-b^2+a+i(2ab+b)\\
&=x_1^2-y_1^2+a+i(2x_1y_1+b)
\end{align}
}
z_2の実部をR、虚部をIとすると
{\begin{align}
z_3 &=(R+Ii)^2+a+ib\\
&=R^2-I^2+a+i(2RI+b)\\
&=x_2^2-y_2^2+a+i(2x_2y_2+b)
\end{align}
}
以上より
{\displaystyle 
\begin{eqnarray}
  \left\{
    \begin{array}{l}
      x_{n+1} = x_n^2-y_2^2+a\\
      y_{n+1} = 2x_ny_n+b
    \end{array}
  \right.
\end{eqnarray}
}
であると推測できる。
証明は分からないので省きます。


プログラム

上記の定義を元にマンデルブロ集合を描いていきます。

仕様
・縦軸を虚数、横軸を実数とする
・範囲は-2~2
・画面中央を原点とする
・縦軸下向きを虚軸の正、横軸右向きを実軸の正とする
・発散の速さに応じて色を変える

収束の速さと対応する色
let num_max = 100; //発散判定のループ回数最大値
let divergence_max = 200; //発散したとみなす基準値

function setup(){
  createCanvas(200, 200);
  colorMode(HSB, 360, 100, 100);
  background(0);
  translate(width/2, height/2); //座標原点を中央に
  scale(width/4, height/4); //座標を-2~2の範囲に
  let xp = 4/width; //ループ1回毎に加算する値(横x軸)
  let yp = 4/height; //ループ1回毎に加算する値(縦y軸)
  noStroke();
  for(let x=-2; x<=2+xp; x+=xp){
    for(let y=-2; y<=2+yp; y+=yp){
      let num = divergence_num(x, y); //発散判定
      if(num==num_max) continue; //発散しなかった場合点を打たない
      num = map(num, 0, num_max, 0, 360, true) //発散の速さに応じて色を変える
      fill(num,100,100);
      rect(x, y, xp, yp);
    }
  }
}

function divergence_num(real, imaginary){ //発散判定
  let n;
  let r = real;
  let i = imaginary;
  for(n=1; n<num_max; n++){
    if(dist(0,0,r,i) > divergence_max) return n; //(Re,Im)の距離が発散した場合nを返す
    let r_tmp = r;
    let i_tmp = i;
    r = pow(r_tmp,2) - pow(i_tmp,2) + real;
    i = 2 * r_tmp * i_tmp + imaginary;
  }
  return n; //発散しなかった場合num_maxを返す
}


プログラムについて

※このプログラムではキャンバスサイズの大きさによって計算量が増えます。
 そのため、サイズを大きくすると細かく描画されますが、高負荷になります。

  • 最初にある変数num_maxとdivergence_maxはそれぞれ数列がどのくらいで発散するか、どのくらいまで計算するかを指定しています。適当に調整してください。
  • map(num, 0, num_max, 0, 360, true)は0~num_maxの範囲を取るnumの値を、0~360の値になるよう変換します。6番目の引数をtrueとすることで、numがnum_maxとなった場合でも360より大きい値を取らなくなります。
  • dist(0,0,r,i)で(0,0)と(r,i)の2点間の距離計算します。複素数の大きさは|a+ib|=\sqrt{a^2+b^2}なので、distで計算することができます。



備考
divergence_num()の計算部分で関数2つの任意の定数CrCiをおき、
{\displaystyle 
\begin{eqnarray}
  \left\{
    \begin{array}{l}
      x_{n+1} = x_n^2-y_2^2+C_r\\
      y_{n+1} = 2x_ny_n+C_i
    \end{array}
  \right.
\end{eqnarray}
}
とすると下図のようなジュリア集合を描画することができます。

ジュリア集合(Cr=-0.034,Ci=-0.710)