2D model plazov
reference: [1] http://guava.physics.uiuc.edu/~nigel/courses/563/Essays_2012/PDF/banerjee.pdf [2] http://blairgemmer.com/docs/BakTang/SandpileWriteup.pdf [3] https://www.youtube.com/watch?v=1MtEUErz7Gg
Izdelali bomo model peščenig plazov v dveh dimenzijah. V ta namen ustvarimo NxN naklonsko matriko naravnih števil (S). Glede na vrednost, ki jo ij-ti element te matrike zavzema, opredelimo S[i,j] kot: -> stabilen S[i,j]<4 -> nestabilen S[i,j]>=4
Na začetku lahko elementi matrike zauzemajo vrednost 0 ali pa njihovo vrednost določimo naključno. Nato pričnemo dodajati zrna peska. Če določena celica postane nestabilna uporabimo naslednjo pravilo: 1) s[i,j]=s[i,j]-4 (celici odvzamemo 4 zrna peska in jih porazdelimo med sosednjimi štirimi celicami) s[i+1,j]=s[i+1,j]+1 s[i-1,j]=s[i-,j]+1 s[i,j+1]=s[i,j+1]+1 s[i,j-1]=s[i,j-1]+1 Ko en element matrike postane nestabilen, uporabimo pravilo (1) in usklajeno z drugimi celicami spremljamo kako velik plaz se sproži (število elementov, ki zaradi pravila (1) posatane nestabilno).
from pylab import* #risanje grafov
from numpy import* #array
from time import sleep #funkcije povezane s casom
import sys #funkcionalnost programa
sys.setrecursionlimit(200000)#stevilo rekurzivnih korakov
import scipy.optimize
def TBW_model(ii,jj,SP):
global event
ip=ii+1
im=ii-1
jp=jj+1
jm=jj-1
if SP[ii,jj]>3:
event+=1
SP[ii,jj]-=4
if ii>0:
SP[im,jj]+=1
TBW_model(im,jj,SP)
if ii<N-1:
SP[ip,jj]+=1
TBW_model(ip,jj,SP)
if jj>0:
SP[ii,jm]+=1
TBW_model(ii,jm,SP)
if jj<N-1:
SP[ii,jp]+=1
TBW_model(ii,jp,SP)
N=21
SP=zeros([N,N],int)
st_korakov=100
plaz=[]
porazdelitev={}
for i in range(st_korakov):
event=0
ii=(N-1)/2
jj=(N-1)/2
SP[ii,jj]+=1
TBW_model(ii,jj,SP)
if event>0:
plaz.append(event)
if event in porazdelitev:
porazdelitev[event]+=1
else:
porazdelitev[event]=1
x=[]
y=[]
for i in porazdelitev:
x.append(i)
y.append(porazdelitev[i])
fig=figure(figsize=(10,10))
ax=fig.add_subplot(111)
ax.scatter(x,y)
ax.set_xscale('log')
ax.set_yscale('log')
def f(x,A,B):
return A*np.power(x,B)
popt, pcov = scipy.optimize.curve_fit(f,x,y)
A, B = popt
perr = np.sqrt(np.diag(pcov))
yy=f(sorted(x),A,B)
p1,=ax.plot(sorted(x),yy,color='r')
ax.legend([p1], [r'${D_f} = %1.3f$'%(B)], loc=1)
show()
