Pythonでカオス・フラクタル
昔，JAVAやらMATLABやらで色々と描いた気がする．ただ，それをそのまま移植すると，ループを回すだけになってしまうので，できる限りPythonicなコードを考えたい．でも，MATLABを使っている時からそうだったけど，ベクトル演算が本当に苦手．色々と参考にしながら，パッと出てくる様になりたい．
参考：
[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])
閃かなかったのでそのまま書き下した．
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')
参考：
[1] Pythonでカオス・フラクタルを見よう！-Qiita
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
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()
広告
ピンバック: フラクタル・ツリー | 粉末@それは風のように (日記)
ピンバック: numpy + PIL | 粉末@それは風のように (日記)