フラクタルの「シダ」

よくフラクタルで登場する「シダ」のグラフィックス。
以下４種のx,yの漸化式が用意されており、確率でどれかを選び、計算を繰り返しながら「点」を描画する。

\begin{equation*}
   \begin{cases}
x_{n+1}=0 \\
y_{n+1}=0.16y_n
   \end{cases}
\end{equation*}

\begin{equation*}
   \begin{cases}
x_{n+1}=0.2x_n-0.26y_n \\
y_{n+1}=0.23x_n+0.22y_n+1.6
   \end{cases}
\end{equation*}

\begin{equation*}
   \begin{cases}
x_{n+1}=-0.15x_n+0.28y_n \\
y_{n+1}=0.26y_n+0.24y_n+0.44
   \end{cases}
\end{equation*}

\begin{equation*}
   \begin{cases}
x_{n+1}=0.85x_n+0.04y_n \\
y_{n+1}=-0.04x_n+0.85y_n+1.6
   \end{cases}
\end{equation*}

回転行列による計算を

\begin{bmatrix}
x_{n+1} \\
y_{n+1} \\
\end{bmatrix}
=
\begin{bmatrix}
cos\theta & -sin\theta \\
sin\theta & cos\theta \\
\end{bmatrix}
\begin{bmatrix}
x_n \\
y_n \\
\end{bmatrix}

とすると

\begin{equation*}
   \begin{cases}
x_{n+1}=x_ncos\theta-y_nsin\theta \\
y_{n+1}=x_nsin\theta+y_ncos\theta
   \end{cases}
\end{equation*}

それで、最後の式で回転角度調べてみようと、

arccos\theta=0.85

\theta=0.554811033rad=0.00968327925degree

0.85^2+0.04^2=0.7241

単純な回転じゃなかった。せん断が入ってるのか？

size(256, 512);
float x=0, y=0;
float x1, y1;
translate(width/2, 0);
for (int n=0; n<100000; n++) {
  float r = random(100);
  if (r<1) {
    x1 = 0;
    y1 = 0.16*y;
  } else if (r<8) {
    x1 =  0.2*x-0.26*y;
    y1 = 0.23*x+0.22*y+1.6;
  } else if (r<15) {
    x1 = -0.15*x+0.28*y;
    y1 =  0.26*x+0.24*y+0.44;
  } else {
    x1 =  0.85*x+0.04*y;
    y1 = -0.04*x+0.85*y+1.6;
  }
  x = x1;
  y = y1;
  stroke(r);
  point(x*40, y*40);
}

シェルピンスキーのギャスケット

0.5倍
0.5倍して(0.25, √3/4)の移動
0.5倍して(0.5, 0)の移動
の3種のx,yの漸化式を用意し、以下同様。

\begin{equation*}
   \begin{cases}
x_{n+1}=0.5x_n \\
y_{n+1}=0.5y_n
   \end{cases}
\end{equation*}

\begin{equation*}
   \begin{cases}
x_{n+1}=0.5x_n+0.25 \\
y_{n+1}=0.5y_n+\frac{\sqrt3}{4}
   \end{cases}
\end{equation*}

\begin{equation*}
   \begin{cases}
x_{n+1}=0.5x_n+0.5 \\
y_{n+1}=0.5y_n
   \end{cases}
\end{equation*}

size(512,512);
float Xn=0, Yn=0;
float X1, Y1;

for (int n=0;n<100000; n++) {
  float r = random(100);
  if (r<33) {
    X1 = 0.5*Xn;
    Y1 = 0.5*Yn;
  } else if (r<66) {
    X1=0.5*Xn+0.25;
    Y1=0.5*Yn+sqrt(3f)/4f;
  } else{
    X1=0.5*Xn+0.5;
    Y1=0.5*Yn;
  }
  Xn=X1;
  Yn=Y1;
  stroke(r);
  point(Xn*512, Yn*512);
}

コッホ曲線

size(512,512);
float Xn=0, Yn=0;
float X1, Y1;
translate(0, height/2);
for (int n=0;n<10000; n++) {
  float r = random(100);
  float s = 1f/3f;
  float t = radians(60);
  if (r<25) {
    X1 = s*Xn;
    Y1 = s*Yn;
  } else if (r<50) {
    X1 = s*cos(t)*Xn - s*sin(t)*Yn + s;
    Y1 = s*sin(t)*Xn + s*cos(t)*Yn;
  } else if (r<75) {
    X1 = s*cos(t)*Xn + s*sin(t)*Yn + cos(t);
    Y1 =-s*sin(t)*Xn + s*cos(t)*Yn + s*sin(t);
  } else{
    X1 = s*Xn + 2*s;
    Y1 = s*Yn;
  }
  Xn = X1;
  Yn = Y1;
  point(Xn*512, Yn*512);
}

ドラゴン曲線

\begin{equation*}
   \begin{cases}
x_{n+1}=0.5x_n-0.5y_n \\
y_{n+1}=0.5x_n+0.5y_n
   \end{cases}
\end{equation*}

\begin{equation*}
   \begin{cases}
x_{n+1}=-0.5x_n-0.5y_n+1 \\
y_{n+1}=0.5x_n-0.5y_n
   \end{cases}
\end{equation*}

size(512,512);
float Xn=0, Yn=0;
float X1, Y1;
translate(width/4, height/2);
for (int n=0;n<200000; n++) {
  float r = random(100);
  if (r<50) {
    X1 = 0.5*Xn-0.5*Yn;
    Y1 = 0.5*Xn+0.5*Yn;
  } else{
    X1 =-0.5*Xn-0.5*Yn + 1;
    Y1 = 0.5*Xn-0.5*Yn;
  }
  Xn = X1;
  Yn = Y1;
  stroke(r);
  point(Xn*256, Yn*256);
}

レヴィドラゴン曲線

\begin{equation*}
   \begin{cases}
x_{n+1}=0.5x_n-0.5y_n \\
y_{n+1}=0.5x_n+0.5y_n
   \end{cases}
\end{equation*}

\begin{equation*}
   \begin{cases}
x_{n+1}=-0.5x_n+0.5y_n+1 \\
y_{n+1}=0.5x_n+0.5y_n
   \end{cases}
\end{equation*}

size(512,512);
float Xn=0, Yn=0;
float X1, Y1;
translate(width/4, height/2);
for (int n=0;n<100000; n++) {
  float r = random(100);
  if (r<50) {
    X1 = 0.5*Xn-0.5*Yn;
    Y1 = 0.5*Xn+0.5*Yn;
  } else{
    X1 =-0.5*Xn+0.5*Yn + 1;
    Y1 = 0.5*Xn+0.5*Yn;
  }
  Xn = X1;
  Yn = Y1;
  stroke(r);
  point(Xn*256, Yn*256);
}

ツインドラゴン

\begin{equation*}
   \begin{cases}
x_{n+1}=0.5x_n-0.5y_n \\
y_{n+1}=0.5x_n+0.5y_n
   \end{cases}
\end{equation*}

\begin{equation*}
   \begin{cases}
x_{n+1}=-0.5x_n+0.5y_n+1 \\
y_{n+1}=-0.5x_n-0.5y_n
   \end{cases}
\end{equation*}

size(512,512);
float Xn=0, Yn=0;
float X1, Y1;
translate(width/4, height/2);
for (int n=0;n<1000000; n++) {
  float r = random(100);
  if (r<50) {
    X1 = 0.5*Xn-0.5*Yn;
    Y1 = 0.5*Xn+0.5*Yn;
  } else{
    X1 =-0.5*Xn+0.5*Yn + 1;
    Y1 =-0.5*Xn-0.5*Yn;
  }
  Xn = X1;
  Yn = Y1;
  stroke(r);
  point(Xn*256, Yn*256);
}

#　参考
反復関数系

Iterated function system