计算物理第十次作业

发表于 更新于 分类于 Waline： 阅读次数：

7.6

Simulate SAWs in three dimensions. Determine the vairation of with step number and find the value of , where this parameter is defined through the relation (7.9). Compare your results with those in Figure 7.6. You should find that decreases for successively higher dimensions. (It is 1 in one dimension and 3/4 in two dimensions.) Can you explain this trend qualitatively?

为了保证随机行走每步的等可能性，这里采用的方法是检查生成的每一个新位置是否与已有位置重复，若重复，则重新生成，直到不重复。为了提高计算精度，模拟1000次后取平均值进行分析。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
import time
import copy
import multiprocessing as mp
import numpy as np
import matplotlib.pyplot as pl
import mpl_toolkits.mplot3d.axes3d as ax3d
from scipy.optimize import curve_fit

def walk(last):
    # 随机挑选三维空间中的一个方向行走
    np.random.seed()
    while True:
        d = np.zeros(3)
        d[np.random.choice([0, 1, 2])] = np.random.choice([-1, 1])
        if (d != last).any():
            return d

def saw(steps = 500):
    while True:
        points = np.zeros((steps + 1, 3))
        points[1] = walk([0,0,0])
        i = 2
        while i <= steps:
            points[i] = points[i - 1] + walk(points[i - 2] - points[i - 1])
            if any(np.equal(points[:i-1], points[i]).all(1)):
                break
            i += 1
        if i == steps + 1:
            r2 = np.power(np.linalg.norm(points, axis=1),2)
            return points, r2
import time
import copy
import multiprocessing as mp
import numpy as np
import matplotlib.pyplot as pl
import mpl_toolkits.mplot3d.axes3d as ax3d
from scipy.optimize import curve_fit

def walk(last):
    # 随机挑选三维空间中的一个方向行走
    np.random.seed()
    while True:
        d = np.zeros(3)
        d[np.random.choice([0, 1, 2])] = np.random.choice([-1, 1])
        if (d != last).any():
            return d

def saw(steps = 50):
    while True:
        points = np.zeros((steps + 1, 3))
        points[1] = walk([0,0,0])
        i = 2
        while i <= steps:
            points[i] = points[i - 1] + walk(points[i - 2] - points[i - 1])
            if any(np.equal(points[:i-1], points[i]).all(1)):
                break
            i += 1
        if i == steps + 1:
            r2 = np.power(np.linalg.norm(points, axis=1),2)
            return points, r2

if __name__ == '__main__':
    t1 = time.time()
    args = [60] * 1000
    pool = mp.Pool()
    results = pool.map(saw, args)
    pool.close()
    pool.join()
    t2 = time.time()
    print('calucte time %.2f s'% (t2-t1))

    r2 = np.array([x[1] for x in results])
    results = np.array([x[0] for x in results])
    r2_mean = np.mean(r2, axis = 0)
    steps = np.arange(len(r2_mean))

    # 对结果进行拟合
    def func(x, a, b):
        return a * np.power(x, b)
    popt, pcov = curve_fit(func, steps, r2_mean)
    fit = func(steps, popt[0], popt[1])

    pl.figure(figsize = (10, 10), dpi = 80)
    pl.plot(steps, r2_mean, 'o', markersize = 1.5, label = '$<r^2>$')
    pl.plot(steps, fit, 'r', label = 'fit')
    pl.text(10, 1, r'$<r^2> \approx %.2f t^{%.3f}$'% (popt[0], popt[1]), size = 20)
    pl.title('Self Avoiding Walk')
    pl.xlabel('Number of Steps')
    pl.ylabel(r'$<r^2>$')
    pl.legend(loc = 'center right')
    pl.show()

    fig = pl.figure(figsize = (10, 10), dpi = 80)
    ax = fig.gca(projection='3d')
    ax.plot(results[0][:,0],results[0][:,1],results[0][:,2])
    ax.plot(results[1][:,0],results[1][:,1],results[1][:,2])
    pl.title('Self Avoiding Walk')
    pl.show()

可以发现拟合的结果是 ，即 ，确实随着维度的增大而减小。直观上解释，可以认为在更高维的情况下可以行走的方向变多，因此在相同的步数下离出发点的距离会变得更短，即路线变得更加弯折。

查找资料可知三维情况下理论值 ，模拟得到的数值和理论值略有偏差，而且重复运行程序会有波动，主要原因是模拟次数和步数不足，使得结果没有足够收敛。但是使用该方法在步数增加的情况下计算量会急剧增大，因为此时绝大多数的模拟情况都会形成重复，可以使用其他动态方法，如爬行算法等，可以大幅减小计算时间。

另外在程序中进行了多进程的尝试，确实可以明显减少计算时间，但要注意应在随机开始时重新生成随机种子，否则会因为进程的复制而使得随机种子相同，造成失去随机的效果。

7.12

Calculate the entropy for the cream-in-your-coffee problem, and reproduce the results in Figure 7.16.

下面采用二维随机行走的方法模拟奶油的扩散，忽略掉粒子的相互作用，并且同一位置可以同时存在多个粒子。对于熵的计算采用单位长度的方格进行统计，即直接统计模拟粒子的坐标分布。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
import copy
import numpy as np
import matplotlib.pyplot as pl

class cream(object):
    def __init__(self, total_time = 200, lenth = 200, particals = 1200):
        self.t = total_time
        self.n = particals
        self.l = lenth
        self.e = []
        self.d = int(self.n / 400)
        max_l = int(self.l / 2 + 10)
        min_l = int(self.l / 2 - 10)
        grid = np.mgrid[min_l:max_l, min_l:max_l]
        points = np.vstack(map(np.ravel, grid)).T
        self.points = np.tile(points, (self.d, 1))

    def run(self):
        t = 0
        diff, counts = np.unique(self.points, axis = 0, return_counts=True)
        self.diff = [copy.deepcopy(diff)]
        self.counts = [copy.deepcopy(counts)]
        while t < self.t:
            old = copy.deepcopy(self.points)
            for i in range(self.n):
                d = np.zeros(2)
                d[np.random.choice((0, 1))] = np.random.choice((-1, 1))
                self.points[i] = old[i] + d
            t += 1
            diff, counts = np.unique(self.points, axis = 0, return_counts=True)
            self.e.append(np.sum(- counts / self.n * np.log(counts / self.n)))
            if t in [100, 500, 1000]:
                self.diff.append(copy.deepcopy(diff))
                self.counts.append(copy.deepcopy(counts))

    def show(self):
        pl.figure(figsize = (10, 10), dpi = 80)
        k = 221
        i = 0
        for t in [0, 100, 500, 1000]:
            pl.subplot(k)
            pl.scatter(self.diff[i][:,0], self.diff[i][:,1], c = self.counts[i], cmap='rainbow', s = 2)
            pl.xlim(0, 200)
            pl.ylim(0, 200)
            pl.title('Cream in coffee t = %d'% t)
            pl.xlabel('x')
            pl.ylabel('y')
            k += 1
            i += 1
        pl.colorbar(orientation='horizontal')
        pl.show()

        pl.figure(figsize = (10, 10), dpi = 80)
        pl.plot(range(self.t), self.e)
        pl.xlabel('time(steps)')
        pl.ylabel('Entropy')
        pl.title('Entropy versus time')
        pl.show()

if __name__ == '__main__':
    a = cream(total_time = 2000, particals=800)
    a.run()
    a.show()

可以看到模拟计算得到的熵值是逐渐增大，并且增长趋势是减缓的，这和课本以及热统的知识相符，即系统会达到熵最大的平衡状态。由于划分的网格比较细，因此可以发现熵的波动幅度大于课本上的图7.16（使用的是8×8方格）。