1

我在更新“reversalsize”列表时遇到问题。它应该每次都更新为“反转”。我需要制作一个图表来检查幂律行为。这是一个简单的元胞自动机，沿通道模拟雪崩行为。我总是在列表中得到相同的值，而不是它应该是一个增加的数字

reversal=0           
reversalsize=[]      
def avalanche(M,KH,pmt):          
p_scale_factor = 5      
global reversal
M_av = M.copy()
M_av_app = M.copy()
for i in range(nx-1):
    if i<nx/3:       
        for j in range(ny-1):
            ix = i+1
            iy = j+1
            p_drop = np.tanh((M[ix,iy]-pmt)/p_scale_factor)   
            rand_drop = np.random.random()
            if rand_drop > 1-p_drop:
                if j==0:    
                    if KH[ix,iy] > pmt:
                        av_mass = (M[ix,iy]-pmt)*ppe                  
                        M_av[ix+1,iy] = M_av_app[ix+1,iy]+3*av_mass/4
                        M_av[ix+1,iy+1] = M_av_app[ix+1,iy+1]+av_mass/4
                        M_av[ix,iy] = M_av_app[ix,iy]-av_mass


                elif 0<j<ny-3:
                    if KH[ix,iy] > pmt:
                        av_mass = (M[ix,iy]-pmt)*ppe
                        M_av[ix+1,iy-1] = M_av_app[ix+1,iy-1]+av_mass/4
                        M_av[ix+1,iy] = M_av_app[ix+1,iy]+av_mass/2
                        M_av[ix+1,iy+1] = M_av_app[ix+1,iy+1]+av_mass/4
                        M_av[ix,iy] = M_av_app[ix,iy]-av_mass

                elif j==ny-3:
                    if KH[ix,iy] > pmt:
                        av_mass = (M[ix,iy]-pmt)*ppe
                        M_av[ix+1,iy] = M_av_app[ix+1,iy]+3*av_mass/4
                        M_av[ix+1,iy-1] = M_av_app[ix+1,iy-1]+av_mass/4
                        M_av[ix,iy] = M_av_app[ix,iy]-av_mass 

                M_av_app = M_av
                reversal +=1
                #print(reversal)
    else:
        ilim1 = i-nx/3        
        ilim2 = nx-(i-nx/3)   
        for j in range(ny-1):
            ix = i+1
            iy = j+1
            p_drop = np.tanh((M[ix,iy]-pmt)/p_scale_factor)
            rand_drop = np.random.random()
            if rand_drop > 1-p_drop:
                if j<ilim1:
                    pass
                if j==ilim1:
                    if KH[ix,iy] > pmt:
                        av_mass = (M[ix,iy]-pmt)*ppe
                        M_av[ix+1,iy+1] = M_av_app[ix+1,iy+1]+av_mass
                        M_av[ix,iy] = M_av_app[ix,iy]-av_mass 

                elif ilim1<j<ilim2:
                    if KH[ix,iy] > pmt:
                        av_mass = (M[ix,iy]-pmt)*ppe
                        M_av[ix+1,iy-1] = M_av_app[ix+1,iy-1]+av_mass/4
                        M_av[ix+1,iy] = M_av_app[ix+1,iy]+av_mass/2
                        M_av[ix+1,iy+1] = M_av_app[ix+1,iy+1]+av_mass/4
                        M_av[ix,iy] = M_av_app[ix,iy]-av_mass
                elif j==ilim2:
                    if KH[ix,iy] > pmt:
                        av_mass = (M[ix,iy]-pmt)*ppe
                        M_av[ix+1,iy-1] = M_av_app[ix+1,iy-1]+av_mass
                        M_av[ix,iy] = M_av_app[ix,iy]-av_mass 
                M_av_app = M_av 
                reversal +=1
                #print(reversal)

return M_av

for i in range(n_iter):        
reversal=0
avalanche(M,KH,pmt)
reversalsize.append(reversal)
print(reversalsize)


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()
4

1 回答 1

1

每次迭代时，都会将反转重新设置为 0。这就是它没有增加的原因：

for i in range(n_iter):
    reversal=0  # Delete this line
于 2020-05-08T17:58:58.900 回答