フラクタル

フラクタル（英: fractal）は、フランスの数学者ブノワ・マンデルブロが導入した幾何学の概念です。 図形の部分と全体が自己相似になっているものなどをいいます。

フラクタルの作成方法は反復計算などを使います。シンプルな計算式で綺麗なフラクタルを作成できます。

1. コッホ曲線(Koch curve)の例

まず、コッホ曲線の例でフラクタルの作成方法を勉強します。

  • ステップ0：方向のある線分を引く。（図左上）
  • ステップ1：線分を3等分し、中央の線分を1辺とする正三角形を線分の左側に描き、下の辺を消す。（図右上）
  • ステップ2：得られた4つの線分に対して同じ操作を繰り返す。（図左下）
  • ステップ3：得られた16の線分に対して同じ操作を繰り返す。（図右下）

コッホ曲線を作成するめに、上記の作業の中の繰り返している部分を考えます。

繰り返している作業：　

  • [1] 与えられる線分を３等分
  • [2] 中央の線分を1辺とする正三角形をこの線分の左側に描き、下の辺を消す。

Step 1 線分の３等分

与えられる線分を以下の配列で表示します。

$$ \left[ \begin{array}{cc}x1 & y1 \\ x2 & y2 \end{array} \right] $$

３等分後の節点は４になりますが、これを以下のように表示します。

$$ \left[ \begin{array}{cc}x1 & y1 \\ x3 & y3 \\ x4 & y4 \\ x2 & y2 \end{array} \right] $$

新しい節点の座標は次のように計算できます。

$$x3=\frac{2}{3} x_1 + \frac{1}{3} x_2, \quad y3=\frac{2}{3} y_1 + \frac{1}{3} y_2 $$$$x4=\frac{1}{3} x_1 + \frac{2}{3} x_2, \quad y4=\frac{1}{3} y_1 + \frac{2}{3} y_2 $$

以下のOctaveのコードを参考してください。

In [5]:
x1=0;y1=0;
x2=1;y2=0;
node_list = [x1,y1;x2,y2]

x3 = 2/3*x1+1/3*x2;
y3 = 2/3*y1+1/3*y2;

x4 = 1/3*x1+2/3*x2;
y4 = 1/3*y1+2/3*y2;

new_node_list = [ x1,y1;  x3,y3;  x4,y4; x2,y2]

node_list =

   0   0
   1   0

new_node_list =

   0.00000   0.00000
   0.33333   0.00000
   0.66667   0.00000
   1.00000   0.00000

上記のコードに作られているnode_listを描画するために、以下の関数を作ります。

In [6]:
function draw_node_list(node_list,plot_option)
  x_list = node_list(:,1);
  y_list = node_list(:,2);
  plot(x_list, y_list, plot_option)
end
In [7]:
draw_node_list(node_list, 'r')
hold on
draw_node_list(new_node_list, 'bo')
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 gnuplot_plot_1a gnuplot_plot_2a

演習１

以下の線分について、３等分を取って、上記の作業に従って新しい点のリストを作ってください。

  • p_list_1 = [0,0; 1,2]
  • p_list_2 = [1,0; -1,1]
In [13]:
p_list_1 = [0,0;1,2];

x1=p_list_1(1,1);y1=p_list_1(1,2);
x2=p_list_1(2,1);y2=p_list_1(2,2);

x3 = 2/3*x1+1/3*x2;
y3 = 2/3*y1+1/3*y2;

x4 = 1/3*x1+2/3*x2;
y4 = 1/3*y1+2/3*y2;

new_node_list = [ x1,y1;  x3,y3;  x4,y4; x2,y2];

%---------------- Simple method to calculate new points. -------------
p1=p_list_1(1,:); p2=p_list_1(2,:);
p3 = 2/3*p1+1/3*p2;
p4 = 1/3*p1+2/3*p2;
new_node_list = [ p1;  p3;  p4; p2];
%-------------------- END ---------------------------------------------


hold on
draw_node_list(p_list_1, 'r-')
hold on
draw_node_list(new_node_list, 'bo')
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1 gnuplot_plot_1a gnuplot_plot_2a

1.2 正三角形の２つ辺を描く

与えられる線分を1辺とする正三角形を考えます。 線分のデータは以下の形を持っています。

$$ \left[ \begin{array}{cc}x1 & y1 \\ x2 & y2 \end{array} \right] $$

指定される正三角形の(x1,y1),(x2,y2)の以外の頂点(x3,y3)を以下のように算出できます。

  • (x2,y2)を(x1,y1)に関して、反時計周り$\theta=\pi/3$回転します。

即ち、

$$ \left( \begin{array}{c} x3\\ y3 \end{array} \right) =\left( \begin{array}{cc} \cos \theta & - \sin \theta \\ \sin \theta & \cos \theta \end{array} \right) \left( \begin{array}{c} x2-x1\\ y2-y1 \end{array} \right) + \left( \begin{array}{c} x1\\ y1 \end{array} \right) $$
In [14]:
x1=0;y1=0;
x2=1;y2=0;
node_list = [x1,y1;x2,y2];

theta=pi/3;
A=[cos(theta), -sin(theta); sin(theta), cos(theta)];

new_point = A*[x2-x1; y2-y1] + [x1;y1];

x3=new_point(1);
y3=new_point(2);
new_node_list = [x1,y1; x3,y3; x2,y2];

hold on;
draw_node_list(node_list,'b');
draw_node_list(new_node_list,'r');
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 gnuplot_plot_1a gnuplot_plot_2a

関数で作業をまとめます

上記の作業をまとめて、関数create_triangleを定義します。

In [37]:
function new_node_list = create_triangle(line)
    x1 = line(1,1); y1 = line(1,2);
    x2 = line(2,1); y2 = line(2,2);

    theta=pi/3;
    A=[cos(theta), -sin(theta); sin(theta), cos(theta)];
    new_point = A*[x2-x1; y2-y1] + [x1;y1];
    x3=new_point(1);
    y3=new_point(2);
    
    new_node_list = [x1,y1; x3,y3; x2,y2];
end

関数 create_triangleを使って、三角形の２辺を作成します。

In [32]:
x1=0;y1=0;
x2=1;y2=0;
line = [x1,y1;x2,y2];
new_node_list = create_triangle(line);

draw_node_list(line,'b');
draw_node_list(new_node_list,'r');
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 gnuplot_plot_1a

演習２

上記の作業に従って、線分から三角形を作ってください。

  • p_list_1 = [0,0;2,2]
  • p_list_2 = [1,1;2,0]

1.3 三等分と三角形作成のコードのまとめ

与えられる線分（node_listで表す）について、三等分と三角形作成の作業を update_line にまとめます。

In [17]:
function new_node_list = update_line (node_list)
    x1 = node_list(1,1); y1 = node_list(1,2);
    x2 = node_list(2,1); y2 = node_list(2,2);

    x3 = 2/3*x1+1/3*x2;
    y3 = 2/3*y1+1/3*y2;

    x4 = 1/3*x1+2/3*x2;
    y4 = 1/3*y1+2/3*y2;
    
    triangle_edges = create_triangle([x3,y3; x4,y4]);
    new_node_list = [ x1,y1;  triangle_edges;  x2,y2 ];
end

演習3

以下の幾つかの線分を update_line で処理して、その結果を確認しなさい。

  • node_list_1 = [0,0; 1,0]
  • node_list_2 = [0,0; 2,2]
In [35]:
%ここにコードを書いてください。

1.4 複数線分の処理

node_listに格納される各点と次ぎの点がなす線分をupdate_lineで処理します。各線分から得られる新しい４つの線分をnew_node_listに入れて、新しい点のリストが作られます。

In [16]:
function new_node_list = update_node_list(node_list)
    new_node_list = [];
    for k = 1 : (size(node_list,1)-1)
       line = node_list(k:k+1,:);
       new_line = update_line(line);
       new_node_list = [new_node_list ; new_line];
    end
end

演習4

以下の線分処理の回数（N）を順番に１，２，３，４，５にしてから結果を確認してください。

In [46]:
node_list = [0,0; 1,0];
N=1;
for k=1:N
    node_list = update_node_list(node_list);
end
draw_node_list(node_list, 'r-')

axis("equal")
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 gnuplot_plot_1a

1.6 コードのまとめ

In [20]:
function new_node_list = create_triangle(line)
    x1 = line(1,1); y1 = line(1,2);
    x2 = line(2,1); y2 = line(2,2);

    theta=pi/3;
    A=[cos(theta), -sin(theta); sin(theta), cos(theta)];
    new_point = A*[x2-x1; y2-y1] + [x1;y1];
    x3=new_point(1);
    y3=new_point(2);
    
    new_node_list = [x1,y1; x3,y3; x2,y2];
    
end 

function new_node_list = update_line (node_list)
    x1 = node_list(1,1); y1 = node_list(1,2);
    x2 = node_list(2,1); y2 = node_list(2,2);

    x3 = 2/3*x1+1/3*x2;
    y3 = 2/3*y1+1/3*y2;

    x4 = 1/3*x1+2/3*x2;
    y4 = 1/3*y1+2/3*y2;
    
    triangle_edges = create_triangle([x3,y3; x4,y4]);
    new_node_list = [ x1,y1;  triangle_edges;  x2,y2 ];
    
end


function new_node_list = update_node_list(node_list)
    new_node_list = [];
    for k = 1 : (size(node_list,1)-1)
       line = node_list(k:k+1,:);
       new_line = update_line(line);
       new_node_list = [new_node_list ; new_line];
    end
end

%------------------　以上は反復計算を定義します。 ------------------

node_list = [0,0;1,0];
N=2;
for k=1:N
    node_list = update_node_list(node_list);
end
draw_node_list(node_list, 'r-')

axis("equal")
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 gnuplot_plot_1a

演習 5

初期の点のリストを正三角形の頂点について、上記のupdate_node_listで処理してみてください。

正三角形の頂点の座標：

$$ (0,0), (1/2,\sqrt{3}/2), (1,0) $$
In [49]:
node_list = [ 0,0;  0.5,sqrt(3)/2;  1,0; 0,0];

??

axis("equal")
axis off

N =  3
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 gnuplot_plot_1a

レポート課題

Minkowski Sausageのフラクタルは以下のような反復計算を使っています。 コッホ曲線の計算コードを参考して、Minkowski Sausageのフラクタルを作ってください。

以下は正方形から作ったMinkowski Sausageの図です。

ヒント

以下の作業の流れに従って、Minkowski Sausageに必要な線分の変形ができます。

Step 1

線分の2つの節点を$p_1=(x1,y1),p_2=(x2,y2)$とします。$p_1$を始点にして、$p_1p_2$線分上$|p1p3|=1/\sqrt{5}$を満たす$p_3$は以下の式で算出できます。

$$ p3 = p1 + \frac{1}{\sqrt{5}} \frac{(p2-p1)}{|p2-p1|} $$
In [32]:
p1=[0,0];
p2=[1,0];
p3 = p1 + 1/sqrt(5)*(p2-p1)/norm(p2-p1)

p3 =

   0.44721   0.00000

Step 2

$p_1$を中心にしてベクトル$p_1p_3$を時計回り$\theta$角度で回転します。ただし、$\theta$は以下ように算出されます。 $$ \theta = \sin^{-1}(\frac{1}{\sqrt{5}}) $$

In [36]:
theta=-asin(1/sqrt(5));
A=[cos(theta), -sin(theta); sin(theta), cos(theta)];

vec = (p3-p1)';
new_point = A*vec + p1';

node_list = [p1; new_point']

hold on
draw_node_list([p1;p2], 'b-o')
draw_node_list(node_list, 'r-o')

axis("equal")

node_list =

   0.00000   0.00000
   0.40000  -0.20000
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0 0.2 0.4 0.6 0.8 1 gnuplot_plot_1a gnuplot_plot_2a

Step 3

線分$p_1p_2$の上の点$p_4$を描く。

  • $p_2$を始点にして、$p_1p_2$線分上$|p_2p_4|=|p_1p_2|/\sqrt{5}$を満たす$p_4$を算出します。
  • $p_2$を中心にして、ベクトル$p_2p_4$を時計回り$\theta$角度で回転します。
In [ ]: