Pythonでカオス・フラクタル

昔，JAVAやらMATLABやらで色々と描いた気がする．ただ，それをそのまま移植すると，ループを回すだけになってしまうので，できる限りPythonicなコードを考えたい．でも，MATLABを使っている時からそうだったけど，ベクトル演算が本当に苦手．色々と参考にしながら，パッと出てくる様になりたい．

・Mandelbrot集合-Wikipedia

参考：

[1] python Mandelbrot マンデルブロ　を　python で描画する-(仮）タイトル 備忘録

[2] Pythonでカオス・フラクタルを見よう！-Qiita

[3] Python + gnuplot でマンデルブロ集合を描く-剥いだ布団とメモランダム.old

画像（2次元データ）として表したり，等高線で表現したり，散布図として頑張ったり，様々な方法があって面白い．Pythonの場合，複素数はそのまま扱えるし，ブロードキャストを巧く扱えるかどうかが肝になってくるんだろう．x, yをnp.newaxisでブロードキャストする方法を最初考えたけど，そういう時にnp.ogird[]というインスタンス（関数じゃない）が便利らしい．「x = y = np.linspace()」みたいになっている時は，これに置き換えられないか考えれば良さそう．

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

x, y = np.ogrid[-2.0:2.0:1000j,-2.0:2.0:1000j]
c = x + 1j*y
z = 0
for i in range(50):
        z = z**2 + c
mask = np.abs(z) < 4

_, ax = plt.subplots(figsize=(8, 6))
ax.imshow(mask.T,extent=[-2, 2, -2, 2])

・エノン写像-Wikipedia

閃かなかったのでそのまま書き下した．

import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

sns.set_style('white')

fxy = lambda a, b, x, y : [1.0 - a * x**2 + y, b * x]
a = 1.4
b = 0.3
xtemp = ytemp = 0
x = [xtemp]
y = [ytemp]
for i in range(10000):
    xtemp, ytemp = fxy(a, b, xtemp, ytemp)
    x.append(xtemp)
    y.append(ytemp)

_, ax = plt.subplots(figsize=(8, 6))
ax.scatter(x, y, color='c')

・ロジスティック写像-Wikipedia

参考：

[1] Pythonでカオス・フラクタルを見よう！-Qiita

[2] Logistic map-BredeWiki

import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

sns.set_style('white')
def logistic(a):
    x = [0.8]
    for i in range(400):
        x.append(a * x[-1] * (1 - x[-1]))
    return x[-100:]

al = [2.0 + i * 0.002 for i in range(1000)]
_, ax = plt.subplots(figsize=(8, 6))
for a in al:
    ax.plot([a]*len(logistic(a)), logistic(a), "c.", markersize=1.7)

・樹木曲線（フラクタル・ツリー）

Turtleを使ってもすごく簡単そう．

参考：

[1] Fractal Tree-Python3 codes

[2] Fractal tree-rosettacode

import math
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

sns.set_style('white')

_, ax = plt.subplots(figsize=(8, 6))
xs = 250
ys = 0
x = [xs]
y = [ys]
def drawTree(x1, y1, angle, depth):
    if depth:
        x2 = x1 + math.cos(math.radians(angle)) * depth * 10.0
        y2 = y1 + math.sin(math.radians(angle)) * depth * 10.0
        x.append(x2)
        y.append(y2)
        drawTree(x2, y2, angle - 15, depth - 1)
        drawTree(x2, y2, angle + 15, depth - 1)
    else:
        ax.plot(x, y)

drawTree(xs, ys, -90, 12)

失敗．

import math
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

sns.set_style('white')

_, ax = plt.subplots(figsize=(8, 6))
xs = 250
ys = 0
def drawTree(x1, y1, angle, depth):
    if depth:
        x2 = x1 + math.cos(math.radians(angle)) * depth * 10.0
        y2 = y1 + math.sin(math.radians(angle)) * depth * 10.0
        ax.plot([x1, x2], [y1, y2], 'c-')
        drawTree(x2, y2, angle - 15, depth - 1)
        drawTree(x2, y2, angle + 15, depth - 1)

drawTree(xs, ys, -90, 12)

参考[1]が凄い綺麗．コードそのままで，実行結果の画像が以下．深さで色を変えるだけでここまで綺麗になるんだね．後，matplotlibは遅いのか，やっている事は一緒の筈だけど，これの方が遥かに速い．（Pillowが高速というのは知っていたけど，ここまで違うとは思わなかった．）

取り敢えず，深さで色を変えてみた．

import math
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import seaborn as sns
%matplotlib inline

sns.set_style('white')

_, ax = plt.subplots(figsize=(8, 6))
xs = 250
ys = 0
def drawTree(x1, y1, angle, depth):
    if depth:
        x2 = x1 + math.cos(math.radians(angle)) * depth * 10.0
        y2 = y1 + math.sin(math.radians(angle)) * depth * 10.0
        ax.plot([x1, x2], [y1, y2], color=cm.hsv(depth/12.0))
        drawTree(x2, y2, angle - 15, depth - 1)
        drawTree(x2, y2, angle + 15, depth - 1)

drawTree(xs, ys, -90, 12)

今ハマっている春色とか．

PAGE_BREAK: PageBreak

import math
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import seaborn as sns
%matplotlib inline

sns.set_style('dark')

_, ax = plt.subplots(figsize=(8, 6))
ax.invert_yaxis()
ax.set_xticklabels([])
ax.set_yticklabels([])
xs = 250
ys = 0
def drawTree(x1, y1, angle, depth):
    if depth:
        x2 = x1 + math.cos(math.radians(angle)) * depth * 10.0
        y2 = y1 + math.sin(math.radians(angle)) * depth * 10.0
        ax.plot([x1, x2], [y1, y2], color=cm.spring(depth/12.0), linewidth=depth)
        drawTree(x2, y2, angle - 15, depth - 1)
        drawTree(x2, y2, angle + 15, depth - 1)

drawTree(xs, ys, -90, 12)

まとめ

import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import seaborn as sns
%matplotlib inline


def mandelbrot():
    sns.set_style('white')

    x, y = np.ogrid[-2.0:2.0:1000j,-2.0:2.0:1000j]
    c = x + 1j*y
    z = 0
    for i in range(50):
        z = z**2 + c
        mask = np.abs(z) < 4

    fig = plt.figure(1, figsize=(8, 6))
    ax = fig.add_subplot(111)
    ax.imshow(mask.T,extent=[-2, 2, -2, 2])

    plt.show()

def henon_map():
    sns.set_style('white')

    def fxy(a, b, x, y): return 1.0 - a * x**2 + y, b * x
    #pep8的に駄目な書き方
    #fxy = lambda a, b, x, y : [1.0 - a * x**2 + y, b * x]
    a = 1.4
    b = 0.3
    xtemp = ytemp = 0
    x = [xtemp]
    y = [ytemp]
    for i in range(10000):
        xtemp, ytemp = fxy(a, b, xtemp, ytemp)
        x.append(xtemp)
        y.append(ytemp)

    fig = plt.figure(2, figsize=(8, 6))
    ax = fig.add_subplot(111)
    ax.scatter(x, y, color='c')

    plt.show()

def logistic_map():
    sns.set_style('white')
    fig = plt.figure(3, figsize=(8, 6))
    ax = fig.add_subplot(111)

    def logistic(a):
        x = [0.8]
        for i in range(400):
            x.append(a * x[-1] * (1 - x[-1]))
        return x[-100:]

    for a in np.linspace(2.0, 4.0, 1000):
        x = logistic(a)
        ax.plot([a]*len(x), x, "c.", markersize=1.7)

    plt.show()

def fractal_tree():
    sns.set_style('dark')
    fig = plt.figure(4, figsize=(8, 6))
    ax = fig.add_subplot(111)
    ax.invert_yaxis()
    ax.set_xticklabels([])
    ax.set_yticklabels([])

    xs = 250
    ys = 0
    def drawTree(x1, y1, angle, depth):
        if depth:
            x2 = x1 + math.cos(math.radians(angle)) * depth * 10.0
            y2 = y1 + math.sin(math.radians(angle)) * depth * 10.0
            ax.plot([x1, x2], [y1, y2], color=cm.spring(depth/12.0), linewidth=depth)
            drawTree(x2, y2, angle - 15, depth - 1)
            drawTree(x2, y2, angle + 15, depth - 1)

    drawTree(xs, ys, -90, 12)

    plt.show()


def main():
    mandelbrot()
    henon_map()
    logistic_map()
    fractal_tree()


if __name__ == '__main__':
        main()