仮想計算機構

IT業界と無縁な派遣社員のブログ

【Python】反復関数系の例


Sierpinskiの三角形

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

プログラム
import numpy as np
import matplotlib.pyplot as plt

def F(i,x,y):
    if i==0:
        return [x/2,y/2]
    elif i==1:
        return [(x+1)/2,y/2]
    else:
        return [x/2,(y+1)/2]

nstep = 30000
x,y = 0,0
X = []
Y = []

for step in range(nstep):
    i = np.random.randint(0,3)
    x,y = F(i,x,y)
    if step > 20:
        X.append(x)
        Y.append(y)

plt.scatter(X,Y,s=0.1,c='black')
plt.show()
実行結果

f:id:riverta1992:20210625222349p:plain

Barnsleyのシダ

バーンズリーのシダ - Wikipedia

プログラム
import numpy as np
import matplotlib.pyplot as plt

def F(r,x,y):
    if r<1:
        return [0,0.16*y]
    elif r<86:
        return [0.85*x+0.04*y,-0.04*x+0.85*y+1.6]
    elif r<93:
        return [0.2*x-0.26*y,0.23*x+0.22*y+1.6]
    else:
        return [-0.15*x+0.28*y,0.26*x+0.24*y+0.44]

nstep = 100000
x,y = 0,0
X = []
Y = []

for step in range(nstep):
    r = np.random.rand()*100
    x,y = F(r,x,y)
    if step > 20:
        X.append(x)
        Y.append(y)

plt.scatter(X,Y,s=0.1,c='green')
plt.show()
実行結果

f:id:riverta1992:20210625221953p:plain

Koch曲線

コッホ曲線 - Wikipedia

プログラム
import numpy as np
import matplotlib.pyplot as plt

def F(r,x,y):
    if r<25:
        return [x/3, y/3]
    elif r<50:
        return [(x-y*np.sqrt(3)+2)/6, (x*np.sqrt(3)+y)/6]
    elif r<75:
        return [(x+y*np.sqrt(3)+3)/6, (-x*np.sqrt(3)+y+np.sqrt(3))/6]
    else:
        return [(x+2)/3, y/3]

nstep = 10000
x,y = 0,0
X = []
Y = []

for step in range(nstep):
    r = np.random.rand()*100
    x,y = F(r,x,y)
    if step > 20:
        X.append(x)
        Y.append(y)

plt.scatter(X,Y,s=0.1,c='black')
plt.show()
実行結果

f:id:riverta1992:20210625223534p:plain

Levy C曲線

Lévy C curve - Wikipedia

プログラム
import numpy as np
import matplotlib.pyplot as plt

def F(i,x,y):
    if i==0:
        return [0.5*x-0.5*y,0.5*x+0.5*y]
    else:
        return [0.5*x+0.5*y+0.5,-0.5*x+0.5*y+0.5]

nstep = 30000
x,y = 0,0
X = []
Y = []

for step in range(nstep):
    i = np.random.randint(0,2)
    x,y = F(i,x,y)
    if step > 20:
        X.append(x)
        Y.append(y)

plt.scatter(X,Y,s=0.1,c='blue')
plt.show()
実行結果

f:id:riverta1992:20210625230014p:plain