フラクタル(I)

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

フラクタルの作成では反復計算などを使っています。フラクタルの特徴の一つは、シンプルな計算式で綺麗なフラクタルを作成することができます。

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

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

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

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

繰り返している作業：　

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

Step 1 線分の３等分

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

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

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

$$ \left[ \begin{array}{cc}x1 & x3 & x4 & x2 \\ y1 & y3 & y4 & 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 $$

例１

与えられる線分の３等分を算出するPythonのコード例

In [1]:
import numpy as np

p1 = np.array([0,0])
p2 = np.array([1,0])

node_list = np.array([p1,p2]).T

p3 = 2/3*p1+1/3*p2;
p4 = 1/3*p1+2/3*p2;

new_node_list = np.array([p1,p3,p4,p2]).T

print(new_node_list)

[[ 0.          0.33333333  0.66666667  1.        ]
 [ 0.          0.          0.          0.        ]]

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

In [3]:
def draw_node_list(node_list,plot_option):
    import numpy as np
    import matplotlib.pyplot as plt
    plt.rcParams['figure.figsize'] = [10, 10]

    plt.axes().set_aspect('equal')
    if node_list.shape[0]:
        x_list = node_list[0,:];
        y_list = node_list[1,:];
        plt.plot(x_list, y_list, plot_option)
    plt.grid()
In [8]:
draw_node_list(node_list, 'ro-')
draw_node_list(new_node_list, 'bo-')

演習１

例１のコードを関数divide_lineに変形して、以下の二つの線分の３等分を取って、新しい節点のリストを描いてください。

  • p_list_1 = np.array([[0,0], [1,2]])
  • p_list_2 = np.array([[1,2], [2,0]])
In [5]:
#ここで新しい節点を作成してください。
def divide_nodes(node_list):
    new_node_list = np.array([])
    # ???
    
    return new_node_list

p_list_1 = np.array([[0,0], [1,2]]).T
p_list_2 = np.array([[1,2], [2,0]]).T

draw_node_list(p_list_1, 'r-')
draw_node_list(p_list_2, 'b-')

new_node_list_1 = divide_nodes(p_list_1)
new_node_list_2 = divide_nodes(p_list_2)

draw_node_list(new_node_list_1, 'ro-')
draw_node_list(new_node_list_2, 'bo-')

1.2 線分から正三角形を作成

与えられる線分をベースにして、その線分を一つの辺とする正三角形の作成を考えます。

線分のデータは以下の形を持っています。

$$ \left[ \begin{array}{cc}x1 & x2 \\ y1 & 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) $$

計算した新しい節点を使って、以下の新しい節点のリストを作成します。

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

In [6]:
p1 = np.array([0,0])
p2 = np.array([1,0])

node_list = np.array([p1,p2]).T

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

p3 = A@(p2-p1) + p1;

new_node_list = np.array([p1,p3,p2]).T
draw_node_list(node_list,'b');
draw_node_list(new_node_list,'r');

関数で作業をまとめます

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

In [7]:
def create_triangle_vertex(p1,p2):

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

    p3 = A@(p2-p1) + p1;    
    return p3

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

In [8]:
p1 = np.array([0,0])
p2 = np.array([1,0])

p3 = create_triangle_vertex(p1,p2);

new_node_list = np.array([p1,p3,p2]).T

draw_node_list(node_list,'bo-');
draw_node_list(new_node_list,'r-');

演習２

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

  • p_list_1 = np.array([[0,0], [1,2]])
  • p_list_2 = np.array([[1,2], [2,0]])

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

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

In [9]:
def divide_nodes(node_list):
    p1 = node_list[:,0]
    p2 = node_list[:,1]
    
    p3 = 2/3*p1+1/3*p2;
    p4 = 1/3*p1+2/3*p2;

    new_node_list = np.array([p1,p3,p4,p2]).T
    return new_node_list

def create_triangle_vertex(p1,p2):

    theta=np.pi/3;
    A=np.array([[np.cos(theta), -np.sin(theta)],[np.sin(theta), np.cos(theta)]]);
    p3 = A@(p2-p1) + p1;
    return p3

def update_two_node_list(node_list):   
    three_node_list = divide_nodes(node_list)
    
    p1 = three_node_list[:,0]
    p2 = three_node_list[:,1]
    p3 = three_node_list[:,2]
    p4 = three_node_list[:,3]

    p5 = create_triangle_vertex(p2,p3) 
    
    return np.array([p1,p2,p5,p3,p4]).T

2点の線分を更新して、新しい節点のリストを作成し、グラフを描きます。

In [10]:
p1 = np.array([0,0])
p2 = np.array([1,0])
node_list = np.array([p1,p2]).T
new_node_list = update_two_node_list(node_list)

draw_node_list(new_node_list,'b-')

演習3

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

  • p_list_1 = np.array([[0,0], [1,2]])
  • p_list_2 = np.array([[1,2], [2,0]])
In [11]:
#ここにコードを書いてください。

1.4 複数節点のリストの処理

以下の関数update_nodes_listはnode_listに格納される各点と次ぎの点がなす線分をupdate_two_node_listで処理します。

各線分から得られる新しい４つの線分の節点をnew_node_listに入れて、新しい点のリストを作成します。

節点の計算では、となりの二つの線分は同じ節点を共有しているので、new_node_listを更新するとき、節点が重複しないように注意してください。

In [12]:
def update_node_list(node_list):
    
    #節点リストの最初の点を新しい節点リストの1番目の点とします。
    new_node_list = [node_list[:,0]]
    
    for k in range(node_list.shape[1]-1):
        two_node_list = node_list[:,k:(k+2)]        
        tmp_node_list = update_two_node_list(two_node_list)
        
        #節点が重複しないように、4点のみを追加します。
        new_node_list.append( tmp_node_list[:,1])
        new_node_list.append( tmp_node_list[:,2])
        new_node_list.append( tmp_node_list[:,3])
        new_node_list.append( tmp_node_list[:,4])
        
    return np.array(new_node_list).T

演習4

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

In [13]:
p1 = np.array([0,0])
p2 = np.array([1,0])
node_list = np.array([p1,p2]).T

N = 1
for k in range(N):
    node_list = update_node_list(node_list)
draw_node_list(node_list,'b-')

1.6 コードのまとめ

In [14]:
#与えられる節点のリストを描画します。
def draw_node_list(node_list,plot_option):
    import numpy as np
    import matplotlib.pyplot as plt
    plt.rcParams['figure.figsize'] = [12, 12]

    plt.axes().set_aspect('equal')
    if node_list.shape[0]:
        x_list = node_list[0,:];
        y_list = node_list[1,:];
        plt.plot(x_list, y_list, plot_option)
    plt.grid()


#与えられる2点のリスト（線分）の３等分を作成します。
def divide_nodes(node_list):
    p1 = node_list[:,0]
    p2 = node_list[:,1]
    
    p3 = 2/3*p1+1/3*p2;
    p4 = 1/3*p1+2/3*p2;

    new_node_list = np.array([p1,p3,p4,p2]).T
    return new_node_list

#与えられる線分によって、正三角形の第3点を作成します。
def create_triangle_vertex(p1,p2):

    theta=np.pi/3;
    A=np.array([[np.cos(theta), -np.sin(theta)],[np.sin(theta), np.cos(theta)]]);
    p3 = A@(p2-p1) + p1;
    return p3

#与えられる線分から4つの線分の節点リストを作成します。
def update_two_node_list(node_list):   
    three_node_list = divide_nodes(node_list)
    
    p1 = three_node_list[:,0]
    p2 = three_node_list[:,1]
    p3 = three_node_list[:,2]
    p4 = three_node_list[:,3]

    p5 = create_triangle_vertex(p2,p3) 
    
    return np.array([p1,p2,p5,p3,p4]).T

def update_node_list(node_list):
    
    #節点リストの最初の点を新しい節点リストの1番目の点とします。
    new_node_list = [node_list[:,0]]
    
    for k in range(node_list.shape[1]-1):
        two_node_list = node_list[:,k:(k+2)]        
        tmp_node_list = update_two_node_list(two_node_list)
        
        #節点が重複しないように、4点のみを追加します。
        new_node_list.append( tmp_node_list[:,1])
        new_node_list.append( tmp_node_list[:,2])
        new_node_list.append( tmp_node_list[:,3])
        new_node_list.append( tmp_node_list[:,4])
        
    return np.array(new_node_list).T



p1 = np.array([0,0])
p2 = np.array([1,0])
node_list = np.array([p1,p2]).T

N = 1
for k in range(N):
    node_list = update_node_list(node_list)
draw_node_list(node_list,'b-')

レポート課題 1

初期の線分のリストを正三角形の三つの辺として、上記のupdate_node_listで処理してみてください。

正三角形の頂点の座標：

$$ (0,0), (1/2,\sqrt{3}/2), (1,0) $$

In [15]:
p1 = np.array([0,0])
p2 = np.array([0.5, np.sqrt(3)/2.0])
p3 = np.array([1,0])

node_list = np.array([p1,p2,p3,p1]).T
draw_node_list(node_list,'b-')

レポート課題2 (オプション)

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 [ ]:
p1=np.array([0,0]);
p2=np.array([1,0]);
p3 = p1 + 1/np.sqrt(5)*(p2-p1)/np.linalg.norm(p2-p1)
p3

Step 2

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

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$角度で回転します。