十大免费ppt模板免费下载网站,搜易网服务内容,丽江北京网站建设,商丘网吧什么时候恢复营业计算物理专题----随机游走实战
Problem 1 Implement the 3D random walk 拟合线 自旋的 拟合函数#xff08;没有数学意义#xff09; 参数#xff1a;0.627,3.336,0.603#xff0c;-3.234
自由程满足在一定范围内的均匀分布以标准自由程为单位长度#xff0c;…计算物理专题----随机游走实战
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 500def mc_experiment():global Lambdaglobal Collisionglobal NewLocation 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)/(i1))**0.5New[i] d[i]/Path#plt.plot(range(Collision),d/Lambda)return Locationfor i in range(Path):l mc_experiment()print(i)if i49:plt.plot(range(Collision),New/Lambda*10,labelpath50)if i99:plt.plot(range(Collision),New/Lambda*5,labelpath100)if i249:plt.plot(range(Collision),New/Lambda*2,labelpath250)if i499:plt.plot(range(Collision),New/Lambda,labelpath500)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 50def mc_experiment():global Lambdaglobal Collisionglobal Newglobal exceedLocation 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,1exceed)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)/(i1))**0.5New[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/(1exceed),labelstr(exceed))f open(./str(j).txt,wb)pickle.dump(New,f)f.close()New np.zeros(Collision)exceed 0.1plt.legend()
plt.title(d/lambda-collision)
plt.pause(0.01)
plt.savefig(1-b.jpg) P1-c.py
import pickleData []
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 dx_data np.array(range(len(Data[0])))[1:]
y_data Data[0][1:]plt.title(curve_fit)
plt.plot(x_data,y_data,r-.,labelraw data)popt,pcov curve_fit(func,x_data,y_data)
plt.plot(x_data,func(x_data,*popt),b--,labelfit 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 pltlamda1 #平均自由程-步长
N1000 #总步数即每次实验走N步t [i for i in range(1,N1)]def drms(m):drms[]#计算均方根距离 for i in range(1,N1,1): #3d-球坐标系利用角参数\thata\phi 描述其移动走N步rnp.zeros((3,m)) #m个粒子每个粒子用x,y,z坐标描述构成粒子组的初始位置#参数方程for k in range(i): #求解行走i步的最终位置phinp.random.uniform(0,2*np.pi,m) #生成m个随机数costhetanp.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 drmsa drms(50)
b drms(500)
c drms(5000)plt.plot(t,a,o,markersize3,marker,label50-paths,colorr)
plt.plot(t,b,o,markersize3,marker*,label500-paths,colorg)
plt.plot(t,c,o,markersize3,markerx,label5000-paths,colorb)
plt.xlabel(Number of collisions)
plt.ylabel(d/lambda)
plt.plot(t,np.sqrt(t),labelSqrt(N),color b)
plt.legend()
plt.show()P-M-2.py
import numpy as np
import matplotlib.pyplot as pltN1000 #总步数即每次实验走N步t [i for i in range(1,N1)]def drms(m,a):drms[]#计算均方根距离 for i in range(1,N1,1): #3d-球坐标系利用角参数\thata\phi 描述其移动走N步rnp.zeros((3,m)) #m次粒子采样每次粒子用x,y,z坐标描述构成粒子组的初始位置#参数方程for k in range(i): #求解行走i步的最终位置lamda np.random.uniform(a,2-a,1)phinp.random.uniform(0,2*np.pi,m) #生成m个随机数costhetanp.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 drmsa 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,markersize3,marker,label0.1-1.9,colorr)
plt.plot(t,b,o,markersize3,marker*,label0.2-1.8,colorg)
plt.plot(t,c,o,markersize3,markerx,label0.3-1.7,colorb)plt.plot(t,d,o,markersize3,markerx,label0.4-1.6,colorr)
plt.plot(t,e,o,markersize3,marker,label0.5-1.5,colorg)
plt.plot(t,f,o,markersize3,marker*,label0.6-1.7,colorb)plt.plot(t,g,o,markersize3,marker*,label0.7-1.3,colorr)
plt.plot(t,h,o,markersize3,markerx,label0.8-1.2,colorg)
plt.plot(t,i,o,markersize3,marker,label0.9-1.1,colorb)plt.xlabel(Number of collisions)
plt.ylabel(d/lambda)
plt.plot(t,np.sqrt(t),labelSqrt(N),color b)
plt.legend()
plt.show()Problem 3 随机游走的正态性校验
P3.py import matplotlib.pyplot as plt
import numpy as np
import timenp.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(Sumi)[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 timenp.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] 1else: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 timenp.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 npnum_walks 100 # 模拟次数
max_steps 100 # 聚合物链长度
step_size 1 # 随机步长Rs [] # 链的端到端距离列表# 多次重复模拟
for i in range(num_walks):positions np.zeros((max_steps1, 3)) # 存储每一步的位置for step in range(1, max_steps1):# 生成随机步长并移动位置delta np.random.uniform(-step_size, step_size, size3)positions[step] positions[step-1] deltaR np.linalg.norm(positions[-1] - positions[0]) # 计算链的端到端距离Rs.append(R)N np.arange(1,max_steps1)
v np.polyfit(np.log(N), np.log(Rs), deg1)[0] # 拟合直线斜率即为 Flory exponentprint(fFlory exponent {v:.3f})##这段代码使用了 NumPy 库来进行向量化计算 ##并通过多次模拟生成了随机游走聚合物链的构型。最后使用最小二乘法拟合直线斜率来估计 Flory exponent 的值。 ## P4 forge.py
import numpy as np
import matplotlib.pyplot as pltnp.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*Times2plt.plot(Times1,D1,lw2)
plt.plot(Times2,D2,lw2)noise1 np.random.uniform(-0.1,0.1,7)
noise2 np.random.uniform(-0.1,0.1,30)D1 noise1
D2 noise2plt.scatter(Times1,D1,s3)
plt.scatter(Times2,D2,s3)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 timenp.random.seed(0)Ne [100,500,1000,3000,10000,20000,50000,100000]
Re []Num 1000for 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)**2e 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),deg1)[0] # 拟合直线斜率即为 Flory exponent
print(v:,v)P4-b.py
import matplotlib.pyplot as plt
import numpy as np
import timenp.random.seed(0)Num 1000Ne [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(temp0.5)[0]temp2 np.where(temp0.5)[0]SUM[j] sum(Choicex[temp1])**2 sum(Choicey[temp2])**2e time.time()print(round(e-s,2),s)Re.append(sum(SUM)/Num)NUM np.arange(1,Num1)
v np.polyfit(2*np.log(np.array(Ne)),np.log(Re),deg1)[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):resFalsefor i in range(len(listx)):if xllistx[i] and yllisty[i]:#成对判断坐标是否已存在resTruereturn res
#判断是否死点
def Die(x,y,listx,listy):x1x10x2x-10y1y-10y2y10ResJudge(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 countx y 0listx[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):#判断死点count1#计数breakelif 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 - 10if Judge(x,y,listx,listy):y10continueelse:listx.append(x)listy.append(y)turtle.setheading(270)turtle.forward(10)elif r 2:#左x - 10if Judge(x,y,listx,listy):x10continueelse:listx.append(x)listy.append(y)turtle.setheading(180)turtle.forward(10)elif r 3:#上y 10if Judge(x,y,listx,listy):y-10continueelse:listx.append(x)listy.append(y)turtle.setheading(90)turtle.forward(10)
#主程序部分
if __name__ __main__:temp aif tempa:turtle.hideturtle()#隐藏画笔Map(16)Draw(16)turtle.done()elif tempb:turtle.tracer(False)#隐藏动画效果for i in range(10,51): #模拟地图规模变化count0#每次变化对死点计数器初始化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) For each of the method,give the N big enough: 2D Sample Random Walk 2D Traditional Random Walk 2D Self Avoiding Random Walk 3,000 is enough (Error:1e-2) 3,000 is enough (Error:1e-2) 50 is enough (Error:1e-2) 其实考虑到自封闭 完全可以将self-avoiding random walk 控制在1e2-1e3上不选1e1下只是不够精确而言。 (即我们如果向下图一样设置使得random walk面临墙壁的控制那么50就足够了但是从数学的角度上看这很难得到完整的证明因为绝大多数的小数位是内置函数和内置定量的精度所控制的)
