尝试用python研究简单的多粒子体系运动情况

这篇文章提供的这种计算方法其实是我前年第一次做的2333，然后我拿这个东西糊过去了学术英语的小论文和计算物理的报告。不是什么特别高深的活，但是写起来很琐碎。

让我们考虑一个各粒子位置和初速度给定的具有n个粒子的体系，在这个体系下，我们只考虑粒子间两两相互作用，于是，在一个二维平面里，这个体系一共会有2n个自由度。假设相互作用只和粒子的位置有关系，而与时间无关，即对于单个粒子，它具有如下大小的相互作用势：

$A_{total,i}=\Sigma_{j}A_{ij}(\vec{r})\\$

现在我们想要知道这个体系在之后的过程中会如何演化，包括体系整体的能量动量会呈现怎样的变化趋势。对于其中单个粒子，它所具有的Hamiltonian为

$H=\dfrac{1}{2}mv^{2}+A_{total,i}(\vec{r})\\$

则由正则方程我们有

$\dfrac{dp_{i}}{dt}=-\dfrac{\partial H}{\partial q_{i}}\\$

这个方程本质上也就是x,y两个分量上的牛顿第二定律

$\ddot x=F_{x},\ddot y=F_{y}\\$

然后我们一共有50个粒子合计100个方程，联立解一下就可以得到最后的结果啦，是不是很简单？

才怪咧

不会有人真的联立几十个微分方程去解的，尤其是每个方程的参量都独立地依赖于初始条件的情况下。

所以为了能让计算机比较取巧地解出这个问题，我们采用分子动力学模拟最基本的方法——使用匀加速运动近似。

将整个演化过程按照 $\Delta t$ 的间隔划分为足够多段，这样每一段运动在这段时间里就都是一个匀加速运动，并且它们都是一阶马尔可夫的——下一段运动只和上一段运动结束时的速度和位置有关系（我们知道，力学系统的终态只取决于位置参量及其对时间的一阶导数）。

接着，为了使计算机能够进一步执行下去运算，我们将这个运动分解在x和y两个方向上，则x和y上的正负表示了方向。最后，我们只要作出各点的坐标即可得到初末态和粒子的轨迹。

我们挑其中一个匀加速运动来分析我们需要做些什么和得到哪些参量：

假设在这段运动开始某个粒子 $i$ 具有速度 $v_{i,x_{0}},v_{i,y_{0}}$ ,它对某个粒子 $j$ 的距离为 $r_{ij}$ ，其分量分别为 $x_{ij},y_{ij}$ ，相互作用为 $F_{ij}$

则粒子 $i$ 在某一方向上受到粒子 $j$ 的作用为

$F_{ij,x}=F_{ij} \cdot \dfrac{x_{ij}}{r_{ij}} \\$

$F_{ij,y}=F_{ij} \cdot \dfrac{y_{ij}}{r_{ij}} \\$

粒子 $i$ 在这个方向上受到的作用力总和为

$F_{ix}=\Sigma_{j}F{ij,x}\\$

$F_{iy}=\Sigma_{j}F{ij,y}\\$

在这个方向上粒子的下一段运动可以写成

$x_{i,0+\Delta t}=x_0+v_{i,x_0}\Delta t+\dfrac{1}{2}F_{ix}\Delta t^2 \\$

$y_{i,0+\Delta t}=y_0+v_{i,y_0}\Delta t+\dfrac{1}{2}F_{iy}\Delta t^2 \\$

这里我直接用力来代表了加速度。并且我们不应该忘记有

$v_{i,x_{0+\Delta t}}=v_{i,x_0}+F_{ix}\Delta t \\$

$v_{i,y_{0+\Delta t}}=v_{i,y_0}+F_{iy}\Delta t \\$

由此，我们只要更新数据，就可以继续进行下一步计算了。当然，你也可以通过修改运动方程前面的系数/权重来使得结果更加精确。

先展示一下一个已经做好的实例，保守力的形式是平方反比的嘛。

看起来似乎这个体系能量并不守恒，毕竟初始状态下是非稳体系。不过事实上坐标轴上的56到53占的比例不到10%，而且从前面看起来整体也是近似守恒的，变化并不十分大。

接下来看一下源代码，毕竟是我自己瞎写的，所以可能效率不高hhh

整个程序只用到了numpy和matplotlib两个模块。

def get_velocity(vx):
    vx=np.random.uniform(-1,1)
    return vx

def get_distance(x0,x1,y0,y1):
    r=pow((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1),0.5)
    return r

def get_angle(x0,x1):
    delta_x=x0-x1
    return delta_x

def get_tan(x,y):
    if y != 0:
        tan = x/y
    else:
        tan=0
    return tan

def get_force(r):
    for i in range(num):
        if r != 0:
            f=1/(r*r)
        else:
            f=0
    return f

上述四个函数用以获得随机的初速度，距离，角度的余弦值和受力等（tan只是因为最开始写顺手了233）

    for l in range(step):
        for m in range(num):
            for n in range(num):
                r[m,n,l]=get_distance(x[m,l],x[n,l],y[m,l],y[n,l])
                f[m,n,l]=get_force(r[m,n,l])
                potenital[m,n,l]=f[m,n,l]*r[m,n,l]
                delta_x[m,n,l]=get_angle(x[m,l],x[n,l])
                delta_y[m,n,l]=get_angle(y[m,l],y[n,l])
                angle_x[m,n,l]=get_tan(delta_x[m,n,l],r[m,n,l])
                angle_y[m,n,l]=get_tan(delta_y[m,n,l],r[m,n,l])
                f_x[m,n,l]=f[m,n,l]*angle_x[m,n,l]
                f_y[m,n,l]=f[m,n,l]*angle_y[m,n,l]
                sumf_x[m,l]+=f_x[m,n,l]
                sumf_y[m,l]+=f_y[m,n,l]
                sump[m,l]+=potenital[m,n,l]
            x[m,l+1]=x[m,l]+0.5*sumf_x[m,l]*t*t+v_x[m,l]*t
            y[m,l+1]=y[m,l]+0.5*sumf_y[m,l]*t*t+v_y[m,l]*t
            v_x[m,l+1]=v_x[m,l]+sumf_x[m,l]*t
            v_y[m,l+1]=v_y[m,l]+sumf_y[m,l]*t
        if l==step-2:
            break

主要的循环过程，利用几个三阶张量存储了中间过程中会出现的针对单个粒子的各个分量。

最后只要用plt等函数作图就可以了。

这里可能会出现的问题在于粒子的“穿模”，在某个过程中极有可能会出现两个粒子互相穿过去的情况，当然，这可以通过将时间步长减到足够小来避免。

全部代码：

import numpy as np
import matplotlib.pyplot  as plt
step=1000
num=9
t=0.006
delta_t=np.zeros((1,step))
x=np.zeros((num,step))
y=np.zeros((num,step))
v_x=np.zeros((num,step))
v_y=np.zeros((num,step))
r=np.zeros((num,num,step))
f=np.zeros((num,num,step))
sumf_x=np.zeros((num,step))
sumf_y=np.zeros((num,step))
delta_x=np.zeros((num,num,step))
delta_y=np.zeros((num,num,step))
angle_x=np.zeros((num,num,step))
angle_y=np.zeros((num,num,step))
potenital=np.zeros((num,num,step))
f_x=np.zeros((num,num,step))
f_y=np.zeros((num,num,step))
p_x=np.zeros((1,step))
p_y=np.zeros((1,step))
p=np.zeros((1,step))
E_v=np.zeros((1,step))
E_p=np.zeros((1,step))
E=np.zeros((1,step))
sump=np.zeros((num,step))




def get_velocity(vx):
    vx=np.random.uniform(-1,1)
    return vx

def get_distance(x0,x1,y0,y1):
    r=pow((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1),0.5)
    return r

def get_angle(x0,x1):
    delta_x=x0-x1
    return delta_x

def get_tan(x,y):
    if y != 0:
        tan = x/y
    else:
        tan=0
    return tan

def get_force(r):
    for i in range(num):
        if r != 0:
            f=1/(r*r)
        else:
            f=0
    return f

if __name__ == "__main__":
    x[:,0]=[0,1,2,0,1,2,0,1,2]
    y[:,0]=[0,0,0,1,1,1,2,2,2]
    plt.scatter(x,y,c='r')
    #plt.title('oringin')
    #plt.show()
    
    for i in range(num):
        v_x[i,0]=get_velocity(v_x[i,0])
        v_y[i,0]=get_velocity(v_y[i,0])



    for l in range(step):
        for m in range(num):
            for n in range(num):
                r[m,n,l]=get_distance(x[m,l],x[n,l],y[m,l],y[n,l])
                f[m,n,l]=get_force(r[m,n,l])
                potenital[m,n,l]=f[m,n,l]*r[m,n,l]
                delta_x[m,n,l]=get_angle(x[m,l],x[n,l])
                delta_y[m,n,l]=get_angle(y[m,l],y[n,l])
                angle_x[m,n,l]=get_tan(delta_x[m,n,l],r[m,n,l])
                angle_y[m,n,l]=get_tan(delta_y[m,n,l],r[m,n,l])
                f_x[m,n,l]=f[m,n,l]*angle_x[m,n,l]
                f_y[m,n,l]=f[m,n,l]*angle_y[m,n,l]
                sumf_x[m,l]+=f_x[m,n,l]
                sumf_y[m,l]+=f_y[m,n,l]
                sump[m,l]+=potenital[m,n,l]
            x[m,l+1]=x[m,l]+0.5*sumf_x[m,l]*t*t+v_x[m,l]*t
            y[m,l+1]=y[m,l]+0.5*sumf_y[m,l]*t*t+v_y[m,l]*t
            v_x[m,l+1]=v_x[m,l]+sumf_x[m,l]*t
            v_y[m,l+1]=v_y[m,l]+sumf_y[m,l]*t
        if l==step-2:
            break
    plt.scatter(x[:,step-1],y[:,step-1],c='b')
    plt.title('position')
    plt.show()
    
    
    
    plt.scatter(x,y,s=2)
    plt.title('trace')
    plt.show()

    

    p_x=np.sum(v_x,axis=0)
    p_y=np.sum(v_y,axis=0)
    p_x=p_x.reshape(1,-1)
    p_y=p_y.reshape(1,-1)
            #print(j,k)
    
    delta_t[0,1]=t
    for c in range(1,step):
        delta_t[0,c]=t+delta_t[0,c-1]

        

    
    for e in range(step):
        for f in range(num):
            E_v[0,e]+=v_x[f,e]**2+v_y[f,e]**2
            E_p[0,e]+=sump[f,e]
    
    
    E=E_v+E_p
    p=p_x+p_y

    plt.plot(delta_t[0],p_x[0])
    plt.title('movementx')
    plt.show()
    plt.plot(delta_t[0],p_y[0])
    plt.title('movementy')
    plt.show()
    plt.plot(delta_t[0],p[0])
    plt.title('movement')
    plt.show()
    plt.plot(delta_t[0],E_v[0])
    plt.title('kinetic')
    plt.show()
    plt.plot(delta_t[0],E_p[0])
    plt.title('potenital')
    plt.show()
    plt.plot(delta_t[0],E[0])
    plt.title('Energy')
    plt.show()

尝试用python研究简单的多粒子体系运动情况

猜你喜欢