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

前回はこちらwoodhero2357.hatenablog.com

致命的なミス

前回書いたプログラムに結構致命的なミスがありました。。。

def diffusion(self):
        start = self.ini_posi()
        x,y = self.randomwalk(start[0],start[1])
        
        flag = True
        while flag:
            x,y = self.randomwalk(x,y)
            flag = self.contact(x,y)        #ここでFalseになっても…
            flag = self.out_range(x,y)　#ここでまたTrueに戻ってしまう…
            
            if abs(x)+abs(y)>15:
                flag = False

この部分です。
粒子が核にくっついたら次の粒子を発射されるようループを終了しないといけない(flagがFalseになる)のですが、
粒子は境界内にいるので次の行で再びflagがTrueに戻り、ループが続いてしまいます。
これによって、次のステップでどの方向に行っても核の隣になってしまい不必要に核が拡大していました。

初期位置の移動

さて、前回は詳しく説明しませんでしたが計算時間短縮のために用いられるテクニックが2つあります。

粒子が核から遠くに行ってしまうとその粒子を諦めて次に行く。

粒子はランダムウォークしているので必ず核にくっついてくれるとは限りません。
あまりにも核から遠く離れてしまうと、再び核に近づいて接触する確率はかなり低くなるでしょう。
これでは計算時間を無駄に消費し続けることになってしまいます。
よって核からある程度離れた粒子は無視することにしましょう。

粒子の初期位置をある程度核の近くにしておき、核が成長するにつれて初期位置も動かす。

場の境界から粒子を発射させたとしても、核に接触するためには必ず核の近くの領域に入らなければなりません。
例えば、境界線が核を中心に持つ半径 $R$ の円だとして、その内側に同じく核が中心の半径 $r(\le R)$ の円があるとします。
境界線から核に到達するには、どうやっても内側の円を通らなければいけませんよね？
なら、内側から発射しても同じじゃないですか？
てか、外側から内側に行くまでの時間無駄じゃないですか？
ということで、初期位置を核の近くするわけです。
もちろん、核が成長していくと初期位置が近くなりすぎてしまうので、ある程度近づいたら初期位置をアップデートします。

今回作ったプログラム

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

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 (x-i)**2 + (y-j)**2 <= 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 np.sqrt(x**2 + y**2)
    
    def diffusion(self,r):
        start = self.ini_posi(r)
        x,y = self.randomwalk(start[0],start[1])
        
        flag = True
        while flag:
            x,y = self.randomwalk(x,y)
            if self.norm(x,y) > 1.5*r:
                break
            flag = self.contact(x,y)

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

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

main()