Final Exercise：Random System

Python physics 2014301020054

---

Abstract

Chapter 7 Random System

In this chapter we consider a class of systems in which randomness plays a central role. These are called random or stochastic systems. Typically, these systems consist of a very large number of "degrees of freedom," which might be associated with particles(e.g.,in a macroscopic sample of a liquid) or perhaps spins(in a piece of a ferromagnet). Randomness can then arise in several different waays. For example,it might be impossible to obtain complete knowledge of the positions and velocities of all the particles. Or, our system may interact with a thermal reservoir in the complicated manner that is best described using probabilies and averages as familiar from statitcal mechanics. Hence, even though the "underlying" physics of the system may be deterministic, our incomplete knowlegde may force us to a statistcal,i.e.,stochastic, description. Foetunately, as we will see, a statistical description is often extremely useful, and indeed, most appropriate in these problems.

---

Background

In Statistical and Thermal Physics, we have done some simulation like "stp_MDApproachToEquilibriumTwoPartitions":

We alse have discussed the behavior of macroscopic systems by appealing to everyday experience and simple observations. So we know that the macroscopic systems are very similar to the following question we discussed.

---

Main Body

Problem7.1

Calculate the diffusion constant analytically for the random-walk simulations in Figures7.3 and 7.4. These simulations were caarried out in one dimension and on a three-dimensional simple-cubic lattice(respectively). Can you generalize your analytic results to a wider class of lattice types such as a two-dimensional triangular lattice, or a three-dimensional face-centered or body-centered cubic lattice?

(7.2)

7.3 Self-Avoiding Walks

Solution

At first We consider a class of systems in which randomness plays a central role.These are called random or stochastic systems.

1.To simplify the problem, we first only consider the one dimension problem.

Firstly,sketch of a random walk in one dimension. The walker began at x=x0=0,and each step is indicated schematically by a dotted arrow.Here the first step happend to be to the right,while the next three steps were to the left.

To plot two figure about $x$ versus step number and <$x^2$> versus step number. The codes are typed as follows. And the figures can be ploted standardly.

Total codes:

 # -*- coding:utf-8 -*- 
import pylab as pl 
import numpy as np 
import math 
import random 
from mpl_toolkits.mplot3d import Axes3D 
import matplotlib.pyplot as plt 
from matplotlib import animation 
class oneparticle(): 
    def __init__(self,x0=0,dt=0.1,T=10,p1=0.5,p2=0.6,x2ave=0,tot=400):
        self.x1=[x0] 
        self.t1=[0] 
        self.x2 = [x0] 
        self.t2 = [0] 
        self.T=T 
        self.dt=dt 
        self.y=[0] 
        self.num=[0] 
        self.num1=[0] 
        self.num2 = [0] 
        self.p1=p1 
        self.p2=p2 
        self.x2ave=[x2ave] 
        self.x=[0] 
        self.y=[0] 
        self.tot=tot 
    def run(self): 
        xnew1=0 
        xnew2=0 
        for j in range (self.tot): 
            x1=0 
            for i in range (100):    #每人走100步 
                if (random.random ()<0.5): 
                    dx=-1 
                else: 
                    dx=1 
                x1=x1+dx 
                self.x.append(x1**2) 
        for k in range (100): 
            xx=0 
            for l in range(self.tot): #取出第k步每人的位移 
                need = k+1+l*100 
                xx=self.x[need]+xx 
            xneed=xx/self.tot    #求出k步位移平方的均值 
            self.x2ave.append(xneed) 
            self.num.append(self.num[-1]+1) 
        ##以下是两个人的位移统计
        for i in range (100): 
            if (random.random() < 0.5): 
                dx = -1 
            else: 
                dx = 1 
            xnew1=xnew1+dx 
            self.x1.append(xnew1) 
            self.num1.append(self.num1[-1]+1) 
            self.y.append(0) 
        for i in range(100): 
            if (random.random() < 0.5): 
                dx = -1 
            else: 
                dx = 1 
            xnew2 = xnew2 + dx 
            self.x2.append(xnew2) 
            self.num2.append(self.num2[-1] + 1) 
        #print(self.num,self.x) 
    def show(self): 
        pl.subplot(122) 
        pl.plot(self.num, self.x2, 'bo', label='$<x^2>$versus time,n=100') 
        pl.xlabel('step number') 
        pl.ylabel('$<x^2>$') 
        pl.title('Random walk in one dimension') 
        pl.xlim(0, 100) 
        pl.ylim(0, 100) 
        pl.axis('equal') 
        pl.legend(loc='upper left', frameon=False) 
        pl.subplot(121) 
        pl.plot(self.num1, self.x1, 'bo', label='walker1') 
        pl.plot(self.num2, self.x2, 'ro', label='walker2') 
        pl.xlabel('step number') 
        pl.ylabel('x') 
        pl.xlim(0,100) 
        pl.ylim(-20.20) 
        pl.title('Random walk in one dimension') 
        #pl.axis('equal') 
        pl.legend(loc='best',frameon=False) 
        pl.show() 
    def show_trajectory(self): 
        fig = plt.figure() 
        ax = plt.axes(title=('Random Walk'), aspect='equal', autoscale_on=False, xlim=(-1.1, 1.1), ylim=(-1.1, 1.1), xlabel=('x'), ylabel=('y')) 
        line = ax.plot([], [], 'b') 
        point = ax.plot([], [], 'ro', markersize=10) 
        images = [] 
        def init(): 
            line = ax.plot([], [], 'b', markersize=8) 
            point = ax.plot([], [], 'ro', markersize=10) 
            return line, point 
        def anmi(i): 
            ax.clear() 
            line = ax.plot(self.x1[0:(1 * i)], self.y[0:(1 * i)], 'b', markersize=8) 
            point = ax.plot(self.x1[(1 * i - 1):(1* i)], self.y[(1 * i - 1):(1 * i)], 'ro', markersize=10) 
            return line, point 
        anmi = animation.FuncAnimation(fig, anmi, init_func=init, frames=1000, interval=1, blit=False) 
        plt.show() 
a=oneparticle() 
a.run() 
a.show() 
a.show_trajectory()

According to the second picture ,we can find that:

$$<x^2>=2Dt$$

$$\sqrt{<x^2>}\sim\sqrt{t}$$

In fact ,this result D can be obtained analytically.Writing the position after n steps, xn ,as a sum of n separate steps gives

$$x_n=\sum_{i=0}^{n}{si}$$

where $si$ is the displacement of the ith step.For this problem, $si=+-1$ with equal probabilities.We can write that

$$x_n^2=\sum_{i=1}^{n}(\sum_{j=1}^{n}s_{i}s_{j})$$

$$s_{i}s_{j}=\pm1,when,i\neq j$$

$$\therefore<x_n^2>=\sum_{i=1}^{n}s_i^2=n$$

since n is equal to time,this is identical to D=1/2

$$\sqrt{<(x_n^2-<x_n^2>)^2>}=\sqrt{<x_n^4-<x_n^2>^2>}$$

$$<x_n^4>=\sum_{i,j,k,l=1}^{n}{<s_{i}s_{j}s_{k}s_{l}>}$$

$$<x_n^4>=\sum_{i=1}^{n}{s_i^4}+3\sum_{i=1}^{n}{[s_i^2\sum_{j\neq i}{s_j^2}]}=3n^2-2n$$

$$\sqrt{<(x_n^2)^2>-<x_n^2>^2}=\sqrt{2n^2-2n}\sim\sqrt{2}n$$

2.Then we consider the step is random too,x versus step is in the range -1 to 1

we just need to change the expression of dx to

dx=random.uniform(-1,1)

And then we plot the follow pictures after several attempts.

we can obtain different results to show the meaning of "random":

3.Then we consider the three-dimensional random walks.

$$r^2\equiv x^2+y^2+z^2$$

Modefied code as follows.

    for j in range(self.tot):
        x = 0
        y = 0
        z = 0
        r = 1
        for i in range(100):  # 每人走100步
            phi = random.uniform(-math.pi, math.pi)
            theta = random.uniform(0, math.pi)
            dx=r*math.sin(theta)*math.cos(phi)
            dy=r*math.sin(theta)*math.sin(phi)
            dz=r*math.cos(theta)
            x = x + dx
            y = y + dy
            z = z + dz
            r_2=x**2+y**2+z**2
            self.r.append(r_2)
    for k in range(100):
        xx = 0
        for l in range(self.tot):  # 取出第k步每人的位移
            need = k + 1 + l * 100
            xx = self.r[need] + xx
        xneed = xx / self.tot  # 求出k步位移平方的均值
        self.r2ave.append(xneed)
        self.num.append(self.num[-1] + 1)

we can also get the animation of one ball's trajectory in three-dimension.

The animation can be obtained by followed code:

from mpl_toolkits.mplot3d import Axes3D  #生成3D图像的包
import matplotlib.pyplot as plt      #用于绘图
from matplotlib import animation #用于制作动画
##此处省略class中前面部分的代码
    def show_trajectory(self):
        fig = plt.figure()
        ax =Axes3D(fig)
        line1=ax.plot([],[],'b:')
        point1=ax.plot([],[],'bo',markersize=10)
        images=[]
        def init():
            line1 = ax.plot([], [], 'b', markersize=8)
            point1 = ax.plot([], [], 'ro', markersize=10)
            return line1, point1
        def anmi(i):
            ax.clear()
            line1 = ax.plot(self.x[0:(1 * i)], self.y[0:(1 * i)], self.z[0:1*i],'b', markersize=8)
            point1 = ax.plot(self.x[(1 * i - 1):(1 * i)], self.y[(1 * i - 1):(1 * i)], self.z[1 * i -1:1 * i],'ro', markersize=10)
            return line1, point1
        anmi = animation.FuncAnimation(fig, anmi, init_func=init, frames=100, interval=1,
                                       blit=False,repeat=False)
        plt.show()

---

If the path followed by our polymer molecule must not be allowed to intersect itself.Only one segment of the polymer molecule is allowed to occupy any particular region of space. A random walk that is subject to this constraint is called a self-avoiding walk,or SAW.

Codes about calculation:

import matplotlib.pyplot as pl
import numpy as np
import math
import random
class particle(object):#def a class"particle"
    diret = 0
    def __init__(self,x,y,state):
        self.x=x
        self.y=y
        self.state=state
    def moveparticle(self):
        if self.diret==0:
            self.x=self.x+0.5
        elif self.diret==1:
            self.x=self.x-0.5
        elif self.diret==2:
            self.y=self.y-0.5
        elif self.diret==3:
            self.y=self.y+0.5
        else :
            print ('error2')
        self.diret=random.randint(0,3)
        return
class screen(object):#def a class"screen"
    lis=[particle(50.0,50.0,0)]
    state=1
    def chargepoint(self,x,y):
        for i in self.lis:
            if i.x==x and i.y==y:
                return 1
        return 0
    def addparticle(self,x=65.0,y=65.0):#add a new particle
        i=self.chargepoint(x,y)
        if self.state==1:
            if i==0:
                self.lis.append(particle(x,y,1))
            else:
                print ('same particle')
                self.state=0
        else:
            print ('full screen')
        for a in self.lis:
            if a.x>66.0 or a.y>66.0 or a.x<34.0 or a.y<34.0:
                self.state=0
        return
    def delparticle(self,x,y):
        for i in range(len(self.lis)):
            if self.lis[i].x==x and self.lis[i]==y:
                self.lis.pop(i)
                self.state=1
                return
        return
    def showscreen(self):
        x=[]
        y=[]
        for i in self.lis:
            x.append(i.x)
            y.append(i.y)
        pl.plot(x,y,'.')
        pl.show()
    def move(self):#random walk of particle
        t1=self.lis.pop()
        while t1.state==1 and t1.x<75.0 and t1.y<75.0 and t1.x>=25 and t1.y>=25:
            i1=self.chargepoint(t1.x+0.5,t1.y)
            i2=self.chargepoint(t1.x-0.5,t1.y)
            i3=self.chargepoint(t1.x,t1.y-0.5)
            i4=self.chargepoint(t1.x,t1.y+0.5)
            s=[i1,i2,i3,i4]
            if i1==0 and i2==0 and i3==0 and i4==0:
                t1.moveparticle()
            elif random.randint(1,100)>85:
                t1.state=0
            elif i1==1 and i2==1 and i3==1 and i4==1:
                t1.state=0
            else:
                while s[t1.diret]==1:
                    t1.diret=random.randint(0,3)
                t1.moveparticle()
            #print 'move a particle'
            #print t1.x,t1.y
        if t1.state==0 :
            self.lis.append(t1)
if __name__=='__main__':
    sc=screen()
    n=0
    while sc.state==1:
        i=random.randint(70,130)
        x=i/2.0
        r=(15**2-(x-50)**2)**0.5#The new particle can appear where r = 15
        k=int(r*2)
        y1=k/2.0
        if random.randint(1,2)==1:
            y=50+y1
        else:
            y=50-y1
        print (x,y)
        sc.addparticle(x,y)
        print ('generate a new particle')
        print (len(sc.lis),sc.state)
        #s=raw_input(':')
        sc.move()
    print (len(sc.lis))
    sc.showscreen()

After this let's consider the cluter growth models.We make some assumptions as follows:

1.we only consider the growth is singler-layer
2.we make a new particle randomly located in the place outside of condensation nuclei(we consider the distance between the particle and the center of the condensation nuclei r>15).Then the particle walk randomly in two dimensions, when tha particle comes into contact with the condensation nuclei, it will stop with a certain probability (adhesion coefficient) and become one part of the condensation nuclei.
3.When the particle walk farther then a certain disstance,we consider the particle wouldn't join the growth of the condensation nuclei.
4.When the condensation meet the border,the activity finished.

粘着系数=1

粘着系数=0.75

粘着系数=0.50

粘着系数=0.15

We can find that the number of particles is become larger with the decrease of the adhesion coefficient.And the shape of condensation nuclei is becoming more and more smooth.

Then we can use the expression to describe the degree of branching

$$D=\frac{lnN}{lnL}$$

$N$ is the number of the particles; $L$ is the boundary dimension.

---

Conclusion

The situations of particle's activity in one,two or three dimensions have been simulated roughly.And the property of the random systems can be analyzed detailedly.

---

Acknowledgement

• Computational Physics Textbook.

• Many PC tools.

• Roomie's help

---

写在结尾的话

计算物理这门课其实一直都令我感兴趣，只是因为这学期比较忙，我没能花太多的时间去钻研，所以之前的作业，包括期末作业，我都交的比较匆忙，有几次可能直到周一晚才交上来，请老师谅解。
然后这门课也让我明白了计算机解决物理问题的有趣，我等寒假有时间了，想慢慢自己去学习，去做长远的实践。
最后谢谢老师的授课啦~

END