1. 二维沙堆模型

算法规则：

• \(L\times L\) 的正方格子上，随机地分布着沙子，每一格上的沙子数不超过 4 。
• 有一个小孩，在 \([L/2, L/2]\) 这个格子上，一粒粒地增加沙子。
• 任何一个格子的沙子数大于等于 4 的时候，就会发生一次沙崩，这个格子上的沙子数 - 4，四个相邻的格子上的沙子数 +1，即沙子流向四个邻居。
• 如果边缘上的沙子发生沙崩，那么一部分沙子流出整个沙盘，例如，如它有 3 个邻居在格子上，则这次沙崩这些邻居的沙子数分别 +1，还有一粒沙子流失。

2. 意义

• 这个沙堆是一个耗散系统，沙子源源不断地从外界流入，又源源不断地流出，系统内的物质并不守恒。
• 如果沙堆的初态确定，那么后续的所有动作就已经完全确定。所以是确定性的系统。
• 随着沙子的累积，小的沙崩是随时会发生的，大的沙崩需要较多的沙子累积，所以较少发生。
• 大的沙崩是很多小的沙崩的叠加。如何从解析上描述自相似性？是否存在标度不变，如何解析描述？
• 如何进一步得到 power law？

3. 代码

以下是一段 python 代码，可以实现小规模的沙堆模型。

import numpy as np
import matplotlib.pyplot as plt

size = 10
nsand = np.zeros([size, size], dtype='uint8')
for i in range(size):
    for j in range(size):
        nsand[i,j] = np.random.randint(0,4)

def avalanch( nsand, size ):
    for i in range( size ):
        for j in range( size ):
            if nsand[i,j] >=4 :
                nsand[i,j] -= 4;
                if i-1 >= 0: nsand[i-1,j] += 1
                if i+1 < size: nsand[i+1,j] += 1
                if j-1 >= 0: nsand[i,j-1] += 1
                if j+1 < size: nsand[i,j+1] += 1
                return 1
    return 0

nstep = 10000
magnitudes = []
for i in range(nstep):
    print("i=",i)
    nsand[(int)(size/2),(int)(size/2)] += 1; magnitude=0
    while avalanch( nsand, size ) > 0 :
        magnitude +=1
        #print("--------------------------------- \n", nsand)
    if magnitude > 0:
        #print("magnitude = ", magnitude)
        magnitudes.append( magnitude )
#print("magnitudes = ", magnitudes)
x = np.bincount(magnitudes)
data = []
for i in range( x.size ): # 去除次数为０的点
    if x[i] > 0: data.append( [ i, x[i] ] )
xdata = []; ydata= [] #
for t in data:
    xdata.append( np.log10(t[0]) ); ydata.append( np.log10(t[1]) )
plt.plot(xdata, ydata, "o")
plt.xlabel(r"$log_{10}$ magnitude")
plt.ylabel(r"$log_{10}$ (number of avalanches)")
plt.savefig("PowerLaw.png")
plt.show()

4. 结果

把沙崩的震级定义为一次大沙崩中包括多少次小沙崩，然后把震级、发声的次数取了 log 以后画在图上，得到下面的图片。

• 曲线非常紧凑，这说明存在一些规律。
• 曲线并不完全是直的，这说明什么？这个弯曲的意义是什么？