1 
nx =10
ny=10
pmt = 7.
SKH_ev = np.zeros((nx,ny))
H_ev = np.zeros((nx,ny))
SKH_ev [1:nx-1,1:ny-1] = 8*np.ones((nx-2,ny-2))
H_ev [1:nx-1,1:ny-1] = np.ones((nx-2,ny-2))
iterations= 100
reversal=0
reversalsize=[]

def avalanche(SKH_ev,H_ev):
    global reversal
    H_av = H_ev.copy()
    H_av_app = H_ev.copy()
    for i in range (1,nx-1):
        for k in range (1,ny-1):
            if SKH_ev[i,k] >= pmt:  #if cell > critical value
                tot = 2.
                H_av[i+1,k] = H_av_app[i+1,k] + tot             #1
                H_av[i+1,k-1] = H_av_app[i+1,k-1] + tot         #2
                H_av[i,k-1] = H_av_app[i,k-1] + tot             #3
                H_av[i-1,k] = H_av_app[i-1,k] + tot             #4
                H_av[i-1,k+1] = H_av_app[i-1,k+1] + tot         #5
                H_av[i,k+1] = H_av_app[i,k+1] + tot             #6
                H_av[i,k] = H_av_app[i,k] -tot                  #0

            H_av_app = H_av
            reversal +=1
    return H_av


AV = avalanche(SKH_ev,H_ev)


for i in range(iterations):   
    avalanche(SKH_ev,H_ev)
    reversalsize.append(reversal)  

#plot    
x, y = np.histogram(reversalsize, 1000)
plt.figure(figsize=(10,7))
plt.clf()
plt.loglog(y[0:-1],x, 'r.')
plt.title("Avalanche Size Distribution", fontsize=14)
plt.xlabel(r"$\log$" + "(Avalanche Size)", fontsize=12)
plt.ylabel(r"$\log$" + "(Frequency)", fontsize=12)
plt.show()

大家好，这是一个简单的沙堆模型代码。细胞是六边形的，细胞的邻域是六边形的。我想检查幂律趋势，但我认为我在更新 reversalsize 数组时遇到了问题。我得到一个始终具有相同值的直方图，而不是通过尾部累积的点获得正确的趋势。

4

0 回答 0