主页 > 软件开发 >

计算物理专题----随机游走实战

软件开发
2025-08-19 00:39:02

计算物理专题----随机游走实战 Problem 1 Implement the 3D random walk

拟合线

自旋的

拟合函数（没有数学意义）

参数：0.627,3.336,0.603，-3.234

自由程满足在一定范围内的均匀分布以标准自由程为单位长度，可得到均匀分布的统计特征方均根距离与平均自由程的比值满足

P1-a.py

import numpy as np from scipy import stats import matplotlib.pyplot as plt # 设置实验参数 Lambda = 1 Collision = 1000 np.random.seed(2) New = np.zeros(Collision) Path = 500 def mc_experiment(): global Lambda global Collision global New Location = np.zeros((Collision,3)) d = np.zeros(Collision) for i in range(1,Collision): theta = np.random.uniform(0,np.pi) phi = np.random.uniform(0,2*np.pi) Location[i] = Location[i-1] + np.array([Lambda*np.sin(theta)*np.cos(phi),\ Lambda*np.sin(theta)*np.sin(phi),\ Lambda*np.cos(theta)]) Dis = np.array([sum(i**2)**0.5 for i in Location]) for i in range(Collision): d[i] = (sum(Dis[:i]**2)/(i+1))**0.5 New[i] += d[i]/Path #plt.plot(range(Collision),d/Lambda) return Location for i in range(Path): l = mc_experiment() print(i) if i==49: plt.plot(range(Collision),New/Lambda*10,label="path=50") if i==99: plt.plot(range(Collision),New/Lambda*5,label="path=100") if i==249: plt.plot(range(Collision),New/Lambda*2,label="path=250") if i==499: plt.plot(range(Collision),New/Lambda,label="path=500") plt.legend() plt.title("<d>/lambda-collision") plt.pause(0.01) plt.savefig("1-a.jpg")

P1-b.py

import numpy as np from scipy import stats import matplotlib.pyplot as plt import pickle # 设置实验参数 exceed = 0.1 Collision = 1000 np.random.seed(2) New = np.zeros(Collision) Path = 50 def mc_experiment(): global Lambda global Collision global New global exceed Location = np.zeros((Collision,3)) d = np.zeros(Collision) for i in range(1,Collision): theta = np.random.uniform(0,np.pi) phi = np.random.uniform(0,2*np.pi) Lambda = np.random.uniform(1-exceed,1+exceed) Location[i] = Location[i-1] + np.array([Lambda*np.sin(theta)*np.cos(phi),\ Lambda*np.sin(theta)*np.sin(phi),\ Lambda*np.cos(theta)]) Dis = np.array([sum(i**2)**0.5 for i in Location]) for i in range(Collision): d[i] = (sum(Dis[:i]**2)/(i+1))**0.5 New[i] += d[i]/Path #plt.plot(range(Collision),d/Lambda) for j in range(6): for i in range(Path): mc_experiment() print(j,":",i) plt.plot(range(Collision),New/(1+exceed),label=str(exceed)) f = open("./"+str(j)+".txt",'wb') pickle.dump(New,f) f.close() New = np.zeros(Collision) exceed += 0.1 plt.legend() plt.title("<d>/lambda-collision") plt.pause(0.01) plt.savefig("1-b.jpg")

P1-c.py

import pickle Data = [] for i in range(6): f = open("./"+str(i)+".txt",'rb') Data.append(pickle.load(f)) import numpy as np from scipy.optimize import curve_fit import matplotlib.pyplot as plt #确定你想要的函数 def func(x,a,b,c,d): return a * np.log(b*x) + c * x**0.5 + d x_data = np.array(range(len(Data[0])))[1:] y_data = Data[0][1:] plt.title("curve_fit") plt.plot(x_data,y_data,"r-.",label="raw data") popt,pcov = curve_fit(func,x_data,y_data) plt.plot(x_data,func(x_data,*popt),"b--",label="fit first") plt.legend() plt.pause(0.01) plt.savefig("1-c") print("popt 1",end=" ") print(popt) print("pcov 1") print(pcov)

P-M-1.py

import numpy as np import matplotlib.pyplot as plt lamda=1 #平均自由程-步长 N=1000 #总步数，即每次实验走N步 t = [i for i in range(1,N+1)] def drms(m): drms=[] #计算均方根距离： for i in range(1,N+1,1): #3d-球坐标系，利用角参数\thata，\phi 描述其移动，走N步 r=np.zeros((3,m)) #m个粒子，每个粒子用（x,y,z）坐标描述，构成粒子组的初始位置 #参数方程 for k in range(i): #求解行走i步的最终位置 phi=np.random.uniform(0,2*np.pi,m) #生成m个随机数 costheta=np.random.uniform(-1,1,m) #生成m个随机数 r[0]=r[0]+lamda*np.sqrt(1-costheta**2)*np.cos(phi) #粒子组的x坐标 r[1]=r[1]+lamda*np.sqrt(1-costheta**2)*np.sin(phi) #粒子组y坐标 r[2]=r[2]+lamda*costheta #粒子组z坐标 d = np.sum(np.reshape(r**2,((r**2).size))) drms.append(np.sqrt(d/m)) #走i次对应的均方根距离 return drms a = drms(50) b = drms(500) c = drms(5000) plt.plot(t,a,'o',markersize='3',marker='+',label='50-paths',color='r') plt.plot(t,b,'o',markersize='3',marker='*',label='500-paths',color='g') plt.plot(t,c,'o',markersize='3',marker='x',label='5000-paths',color='b') plt.xlabel('Number of collisions') plt.ylabel('<d>/lambda') plt.plot(t,np.sqrt(t),label='Sqrt(N)',color = 'b') plt.legend() plt.show()

P-M-2.py

import numpy as np import matplotlib.pyplot as plt N=1000 #总步数，即每次实验走N步 t = [i for i in range(1,N+1)] def drms(m,a): drms=[] #计算均方根距离： for i in range(1,N+1,1): #3d-球坐标系，利用角参数\thata，\phi 描述其移动，走N步 r=np.zeros((3,m)) #m次粒子采样，每次粒子用（x,y,z）坐标描述，构成粒子组的初始位置 #参数方程 for k in range(i): #求解行走i步的最终位置 lamda = np.random.uniform(a,2-a,1) phi=np.random.uniform(0,2*np.pi,m) #生成m个随机数 costheta=np.random.uniform(-1,1,m) #生成m个随机数 r[0]=r[0]+lamda*np.sqrt(1-costheta**2)*np.cos(phi) #粒子组的x坐标 r[1]=r[1]+lamda*np.sqrt(1-costheta**2)*np.sin(phi) #粒子组y坐标 r[2]=r[2]+lamda*costheta #粒子组z坐标 d = np.sum(np.reshape(r**2,((r**2).size))) drms.append(np.sqrt(d/m)) return drms a = drms(500,0.1) b = drms(500,0.2) c = drms(500,0.3) d = drms(500,0.4) e = drms(500,0.5) f = drms(500,0.6) g = drms(500,0.7) h = drms(500,0.8) i = drms(500,0.9) plt.plot(t,a,'o',markersize='3',marker='+',label='0.1-1.9',color='r') plt.plot(t,b,'o',markersize='3',marker='*',label='0.2-1.8',color='g') plt.plot(t,c,'o',markersize='3',marker='x',label='0.3-1.7',color='b') plt.plot(t,d,'o',markersize='3',marker='x',label='0.4-1.6',color='r') plt.plot(t,e,'o',markersize='3',marker='+',label='0.5-1.5',color='g') plt.plot(t,f,'o',markersize='3',marker='*',label='0.6-1.7',color='b') plt.plot(t,g,'o',markersize='3',marker='*',label='0.7-1.3',color='r') plt.plot(t,h,'o',markersize='3',marker='x',label='0.8-1.2',color='g') plt.plot(t,i,'o',markersize='3',marker='+',label='0.9-1.1',color='b') plt.xlabel('Number of collisions') plt.ylabel('<d>/lambda') plt.plot(t,np.sqrt(t),label='Sqrt(N)',color = 'b') plt.legend() plt.show() Problem 3 随机游走的正态性校验

P3.py

import matplotlib.pyplot as plt import numpy as np import time np.random.seed(0) s = time.time() N = 100000 N = int(N) Num = 10000 Num = int(Num) Choice = np.random.choice([-1,1],(N,Num)) Sum = sum(Choice[:,]) e = time.time() print("time:",round(e-s,2)) ##plt.hist(Sum,50) ##plt.title("Distribution of position") ##plt.savefig("Distribution of position.jpg") ##plt.pause(0.01) Position = np.zeros(2061) for i in range(-1030,1031): Position[i] = len(np.where(Sum>i)[0])/Num ##plt.plot(range(1031),Position) ##plt.savefig("P3-c.jpg") ##plt.pause(0.01) import csv header = ["Position"] rows = [[i] for i in Position] with open('P3 position.csv','w',newline="") as file: writer = csv.writer(file) writer.writerow(header) writer.writerows(rows)

从前面的图中可以看出，对于足够大的N，计算出的分布可以用高斯分布来近似

样本量

中位数

平均值

标准差

偏度

峰度

S-W检验

K-S检验

2061

0.502

0.5

0.405

-0.001

-1.713

0.829(0.000***)

0.149(1.1e-40)

P3-e.py

import matplotlib.pyplot as plt import numpy as np import time np.random.seed(0) s = time.time() #step:N N = 3000 N = int(N) #repeat:Num Num = 10000 Num = int(Num) Choice = np.random.random((N,Num)) CHOICE = np.zeros((N,Num)) for i in range(N): for j in range(Num): if Choice[i][j] <= 0.7: CHOICE[i][j] = 1 else: CHOICE[i][j] = -1 Sum = sum(CHOICE[:,]) e = time.time() print("time:",round(e-s,2)) plt.hist(Sum,50) plt.title("Distribution of position-e") plt.savefig("Distribution of position-e N3000.jpg") plt.pause(0.01) import csv header = ["Position"] rows = [[i] for i in Sum] with open('P3-e position N3000.csv','w',newline="") as file: writer = csv.writer(file) writer.writerow(header) writer.writerows(rows)

修改概率使得向正向移动概率为0.7

P3-f.py

import matplotlib.pyplot as plt import numpy as np import time np.random.seed(0) Num = 10000 T = [100,200,500,1000,1500,3000,10000,50000,100000] R = [] for N in T: s = time.time() Choice = np.random.choice([-1,1],(N,Num)) Sum = sum(Choice[:,]) R.append(sum(Sum**2)/Num) e = time.time() print("time:",round(e-s,2)) plt.loglog(T,R) plt.title("log-log E(x^2)-Num") plt.savefig("P3-f-2.jpg") plt.pause(0.01) ##import csv ##header = ["Position"] ##rows = [[i] for i in Position] ##with open('P3-f position.csv','w',newline="") as file: ## writer = csv.writer(file) ## writer.writerow(header) ## writer.writerows(rows)

走N步,轴上移动的距离为X

Problem 4 二维随机游走的自封闭性

Flory exponent.py

##Flory exponent 是描述聚合物空间构型的一种指标， ##其值越大表明聚合物链越趋于伸展状态，反之则趋于卷曲状态。 ## ##在随机游走模型中， ##可以通过生成随机步长并多次重复步骤来模拟聚合物链的构型演化。 ##通过计算链的端到端距离 $R$ 与聚合物链长度 $N$ 之间的关系，可以得到 Flory exponent $v$ 的估计值。 ##

import numpy as np num_walks = 100 # 模拟次数 max_steps = 100 # 聚合物链长度 step_size = 1 # 随机步长 Rs = [] # 链的端到端距离列表 # 多次重复模拟 for i in range(num_walks): positions = np.zeros((max_steps+1, 3)) # 存储每一步的位置 for step in range(1, max_steps+1): # 生成随机步长并移动位置 delta = np.random.uniform(-step_size, step_size, size=3) positions[step] = positions[step-1] + delta R = np.linalg.norm(positions[-1] - positions[0]) # 计算链的端到端距离 Rs.append(R) N = np.arange(1,max_steps+1) v = np.polyfit(np.log(N), np.log(Rs), deg=1)[0] # 拟合直线斜率即为 Flory exponent print(f"Flory exponent = {v:.3f}")

##这段代码使用了 NumPy 库来进行向量化计算， ##并通过多次模拟生成了随机游走聚合物链的构型。最后，使用最小二乘法拟合直线斜率来估计 Flory exponent 的值。 ##

P4 forge.py

import numpy as np import matplotlib.pyplot as plt np.random.seed(0) Times1 = np.array([0.8,1.1,1.5,1.8,2.0,2.1,2.4]) Times2 = np.linspace(2.5,6,30) D1 = 4/3*Times1 D2 = 4/3*Times2 plt.plot(Times1,D1,lw=2) plt.plot(Times2,D2,lw=2) noise1 = np.random.uniform(-0.1,0.1,7) noise2 = np.random.uniform(-0.1,0.1,30) D1 += noise1 D2 += noise2 plt.scatter(Times1,D1,s=3) plt.scatter(Times2,D2,s=3) plt.xlabel("Time") plt.ylabel("$D^2$") plt.title("<D^2> versus T for self avoiding walk in 2D") plt.pause(0.01)

P4-a.py

import matplotlib.pyplot as plt import numpy as np import time np.random.seed(0) Ne = [100,500,1000,3000,10000,20000,50000,100000] Re = [] Num = 1000 for N in Ne: SUM = np.zeros(Num) s = time.time() for j in range(Num): Choicex = np.random.choice([-1,1],N) Choicey = np.random.choice([-1,1],N) SUM[j] = sum(Choicex)**2 + sum(Choicey)**2 e = time.time() print(round(e-s,2),"s") Re.append(sum(SUM)/Num) ##plt.hist(SUM,50) ##plt.title("Distribution of position 2D sample") ##plt.pause(0.01) v = np.polyfit(2*np.log(np.array(Ne)),np.log(Re),deg=1)[0] # 拟合直线斜率即为 Flory exponent print("v:",v)

P4-b.py

import matplotlib.pyplot as plt import numpy as np import time np.random.seed(0) Num = 1000 Ne = [100,500,1000,3000,10000,20000,50000,100000] Re = [] for N in Ne: SUM = np.zeros(Num) s = time.time() for j in range(Num): Choicex = np.random.choice([-1,1],N) Choicey = np.random.choice([-1,1],N) temp = np.random.random(N) temp1 = np.where(temp>=0.5)[0] temp2 = np.where(temp<0.5)[0] SUM[j] = sum(Choicex[temp1])**2 + sum(Choicey[temp2])**2 e = time.time() print(round(e-s,2),"s") Re.append(sum(SUM)/Num) NUM = np.arange(1,Num+1) v = np.polyfit(2*np.log(np.array(Ne)),np.log(Re),deg=1)[0] # 拟合直线斜率即为 Flory exponent print("v:",v) ##plt.hist(SUM,50) ##plt.title("Distribution of position 2D sample") ##plt.pause(0.01)

P4-图像绘制.py

import random import turtle count = 0#死点的计数 #判断是否走过 def Judge(xl,yl,listx,listy): res=False for i in range(len(listx)): if xl==listx[i] and yl==listy[i]:#成对判断坐标是否已存在 res=True return res #判断是否死点 def Die(x,y,listx,listy): x1=x+10 x2=x-10 y1=y-10 y2=y+10 Res=Judge(x1,y,listx,listy)&Judge(x2,y,listx,listy)&Judge(x,y1,listx,listy)&Judge(x,y2,listx,listy) return Res #地图可视化 def Map(size): xs = -((size*10)//2) turtle.pensize(1) turtle.speed(10) #纵坐标的线绘制 for y in range(-((size*10)//2),((size*10)//2)+1,10): turtle.penup() turtle.goto(xs,y) turtle.pendown() turtle.forward(size*10) #横坐标线绘制 ys = ((size*10)//2) turtle.right(90) for x in range(-((size*10)//2),((size*10)//2)+1,10): turtle.penup() turtle.goto(x,ys) turtle.pendown() turtle.forward(size*10) #路径绘制函数 def Draw(size): global count x = y = 0 listx=[0] listy=[0] #设定笔的属性 turtle.pensize(2) turtle.speed(0) turtle.color("red") #模拟走动(是个方向等概率） turtle.penup() turtle.goto(0,0) turtle.pendown() while abs(x) < ((size*10)//2) and abs(y) < ((size*10)//2): r = random.randint(0,3)#产生随机数，0右，1下，2左，3上表示是个方向 if Die(x,y,listx,listy):#判断死点 count+=1#计数 break elif r == 0:#右 x += 10 if Judge(x,y,listx,listy):#判断是否为走过的点 x-=10 #是的话坐标不变 continue#终止本次循环 else: listx.append(x) listy.append(y) turtle.setheading(0) turtle.forward(10) elif r == 1:#下 y -= 10 if Judge(x,y,listx,listy): y+=10 continue else: listx.append(x) listy.append(y) turtle.setheading(270) turtle.forward(10) elif r == 2:#左 x -= 10 if Judge(x,y,listx,listy): x+=10 continue else: listx.append(x) listy.append(y) turtle.setheading(180) turtle.forward(10) elif r == 3:#上 y += 10 if Judge(x,y,listx,listy): y-=10 continue else: listx.append(x) listy.append(y) turtle.setheading(90) turtle.forward(10) #主程序部分 if __name__ == "__main__": temp = 'a' if temp=='a': turtle.hideturtle()#隐藏画笔 Map(16) Draw(16) turtle.done() elif temp=='b': turtle.tracer(False)#隐藏动画效果 for i in range(10,51): #模拟地图规模变化 count=0#每次变化对死点计数器初始化 for j in range(0,10000):#10000次仿真训练 Draw(i) turtle.reset() print('For lattice of size ',i,', the probability of dead-end paths is ',count/100,'%') else: print('input error')

2D Sample Random Walk

拟合直线斜率

v: 0.5022164965587219

选取点

100,500,1000,3000,10000,20000,50000,100000

2D Traditional Random Walk

选取点 100,500,1000,3000,10000,20000,50000,100000

拟合直线斜率 v: 0.49883658055370034

2D Self-Avoiding Random Walk

选取点 Range(2,20)

拟合直线1斜率 v: 1.3074916500876987

拟合直线2斜率 v: 1.502393127(3/4*2)