PythonでDLAを描く(3) - マンデルブログ

まとめ

けっこう飽きてきました(笑)
細かい修正を加えて、終わりにしようと思います。

import random
import math
import numpy as np
import matplotlib.pyplot as plt
import time

class Diffusion():
    def __init__(self):
        self.seed = [[0,0]]
        
    def randomwalk(self,x,y):
        a = random.choice([[1,0],[-1,0],[0,1],[0,-1]])
        return x+a[0],y+a[1]
    
    def contact(self,x,y):
        for i,j in self.seed:
            if self.norm(x-i,y-j) == 1:
                self.seed.append([x,y])
                return False
        return True
    
    def ini_posi(self,r):
        theta = 2 * np.pi * np.random.rand()
        i = int(r*np.cos(theta))
        j = int(r*np.sin(theta))
        return i,j
    
    def norm(self,x,y):
        return math.sqrt(x**2+y**2)
    
    def reset(self,r):
        return self.ini_posi(r)
    
    def diffusion(self,r):
        ini_x,ini_y = self.ini_posi(r)
        x,y = self.randomwalk(ini_x,ini_y)
        
        flag = True
        while flag:
            x,y = self.randomwalk(x,y)
            if self.norm(x,y) > 1.5*r:
                x,y = self.ini_posi(r)
            flag = self.contact(x,y)

def main():
    d = Diffusion()
    r = 5
    for n in range(500):
        d.diffusion(r)
        if n % 3*r == 0:
            norm_list = [d.norm(i,j) for i,j in d.seed]
            if np.max(norm_list)>= r-1:
                r = 2*r

    print(len(d.seed))
    plt.figure(figsize=(8,8))
    plt.plot(*np.array(d.seed).T,".")
    plt.savefig("DLA.png", bbox_inches="tight", pad_inches=0.0)

start = time.time()
main()
elapsed_time = time.time() - start
print("elapsed time:",elapsed_time,"s")

主な変更点

DLAでは粒子数 $N$ に対して、 $\displaystyle N\propto R^{d_D}$ が成り立ちます。 $R$ がDLAクラスターの回転半径、 $d_D$ がフラクタル次元です。DLAのフラクタル次元はだいたい1.7くらいなので、粒子が回転半径のだいたい3倍くらいになったら境界をチェックすることにしました。