1. H=p22+x22H=\frac{p^2}{2}+\frac{x^2}{2}

    可知,在不考虑噪声、不考虑粘滞阻力的情况下,粒子的运动规律为:

    x¨=p˙=Hx=x\ddot{x}=\dot{p}=-\frac{\partial H}{\partial x}=-x

    对上式加入噪声和粘滞,即得此时的郎之万方程:

    x¨+γx˙+x=Γ(t)\ddot{x}+\gamma\dot{x}+x=\Gamma(t)

  2. 对于白噪声,有:

    Γ(t)=0Γ(t)Γ(t)=2Dδ(tt)\begin{aligned} \langle\Gamma(t)\rangle&=0 \\ \langle\Gamma(t)\Gamma(t')\rangle&=2D\delta(t-t') \end{aligned}

    其中,DD 满足涨落耗散定理:

    D=γkBTmD=\dfrac{\gamma k_BT}{m}

    加上郎之万方程:

    x¨+γx˙f(x)=Γ(t)\ddot{x}+\gamma\dot{x}-f(x)=\Gamma(t)

    就可以推导得到态密度函数:

    ρ(p,x)=NeH(p,x)/(kBT)\rho(p, x)=Ne^{H(p, x)/(k_BT)}

    其中,NN 为归一化常数。

    对于这道题,取 kB=1k_B=1,也就是:

    ρ(p,x)=12πTe(x2+p2)/(2T)\rho(p, x)=\frac{1}{2\pi T}e^{-(x^2+p^2)/(2T)}

    那么,就可以推算:

    E=H=12x2+p2=12+(x2+p2)ρ(p,x)dxdp=14πT+(x2+p2)e(x2+p2)/(2T)dxdp=14πT(+x2ex2/(2T)dx+ep2/(2T)dp++p2ep2/(2T)dp+ex2/(2T)dx)=12πT+x2ex2/(2T)dx=T\begin{aligned} \overline{E}&=\langle H \rangle \\ &=\dfrac{1}{2}\langle x^2+p^2\rangle \\ &=\dfrac{1}{2}\iint_{-\infty}^{+\infty}(x^2+p^2)\rho(p, x)\text{d}x\text{d}p \\ &=\dfrac{1}{4\pi T}\iint_{-\infty}^{+\infty}(x^2+p^2)e^{-(x^2+p^2)/(2T)}\text{d}x\text{d}p \\ &=\dfrac{1}{4\pi T}\left(\int_{-\infty}^{+\infty}x^2e^{-x^2/(2T)}\text{d}x\int_{-\infty}^{+\infty}e^{-p^2/(2T)}\text{d}p+\int_{-\infty}^{+\infty}p^2e^{-p^2/(2T)}\text{d}p\int_{-\infty}^{+\infty}e^{-x^2/(2T)}\text{d}x\right) \\ &=\dfrac{1}{\sqrt{2\pi T}}\int_{-\infty}^{+\infty}x^2e^{-x^2/(2T)}\text{d}x \\ &=T \end{aligned}

  3. 程序的思路是这样的:

    1. 生成一堆粒子,具有随机的初始位置和动量。
    2. 给定不同的温度 TT,进行较长时间的演化。
    3. 计算末态的能量,与第二问结果对比。

    演化满足郎之万方程:

    x¨+γx˙+x=Γ(t)\ddot{x}+\gamma\dot{x}+x=\Gamma(t)

    也即:

    x˙=pp˙=γpxΓ(t)\begin{aligned} \dot{x}&=p \\ \dot{p}&=-\gamma p-x-\Gamma(t) \end{aligned}

    假定每步演化的时间为 hh。一阶近似的结果是:

    xnewxold=fx(xold,pold)hpnewpold=fy(xold,pold)h+Γ\begin{aligned} x_\text{new}-x_\text{old}&=f_x(x_\text{old}, p_\text{old})h \\ p_\text{new}-p_\text{old}&=f_y(x_\text{old}, p_\text{old})h + \Gamma \end{aligned}

    其中,fxf_xfpf_p 的定义为:

    fx(x,p)=pfp(x,p)=γpx\begin{aligned} f_x(x, p)&=p \\ f_p(x, p)&=-\gamma p-x \end{aligned}

    Γ\Gamma 的定义为:

    Γ=2γThξ\Gamma=\sqrt{2\gamma Th}\xi

    其中,ξ\xi 为每次迭代时生成的一个高斯随机数。γ\gamma 为粘滞系数。

    我们使用二阶近似(二阶龙格库塔算法):

    xnewxold=12(fx(xold,pold)+fx(xnew,pnew))hpnewpold=12(fp(xold,pold)+fp(xnew,pnew))h+Γ\begin{aligned} x_\text{new}-x_\text{old}&=\frac{1}{2}(f_x(x_\text{old}, p_\text{old})+f_x(x'_\text{new}, p'_\text{new}))h \\ p_\text{new}-p_\text{old}&=\dfrac{1}{2}(f_p(x_\text{old}, p_\text{old})+f_p(x'_\text{new}, p'_\text{new}))h + \Gamma \end{aligned}

    其中,xnewx'_\text{new}pnewp'_\text{new} 的定义为:

    xnew=xold+fx(xold,pold)hpnew=pold+fp(xold,pold)h+Γ\begin{aligned} x'_\text{new}&=x_\text{old}+f_x(x_\text{old},p_\text{old})h \\ p'_\text{new}&=p_\text{old}+f_p(x_\text{old},p_\text{old})h+\Gamma \end{aligned}

    在同一次迭代中,Γ\Gamma 取相同的值。

    取参数:

    n=104, h=102, t=100, γ=0.1, T[0,100]n=10^4,\ h=10^{-2},\ t=100,\ \gamma=0.1,\ T\in[0, 100]

    得到了以下结果:

    Figure_1

    可以看到,在温度较低的时候,实验与模拟符合得很好;在温度较高时,会出现一些偏差。对于高温区,放大来看,是这样的:

    Figure_4

    再增大粒子数去模拟。取参数:

    n=105, h=102, t=100, γ=0.1, T[90,100]n=10^5,\ h=10^{-2},\ t=100,\ \gamma=0.1,\ T\in[90, 100]

    Figure_2

    减小步长:

    n=104, h=103, t=100, γ=0.1, T[90,100]n=10^4,\ h=10^{-3},\ t=100,\ \gamma=0.1,\ T\in[90, 100]

    Figure_3

    增加模拟时间:

    n=104, h=102, t=103, γ=0.1, T[90,100]n=10^4,\ h=10^{-2},\ t=10^3,\ \gamma=0.1,\ T\in[90, 100]

    Figure_5

    可以看到:当温度增加时,确实会使得系统更加不稳定;而这可以通过增加粒子数(系综数)来解决;增加模拟时间和减小步长用处不大。

    代码:

    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
    import numpy
    import matplotlib.pyplot

    numpy.random.seed(123)
    n = int(1e4) # 粒子数
    h = 1e-2 # 步长
    t_end = 1000 # 演化时间
    gamma = 0.1 # 粘滞系数
    T = numpy.linspace(90, 100, 11).reshape([11, 1]) # 模拟的温度范围
    x = numpy.random.rand(n, 1)
    p = numpy.random.rand(n, 1)
    H = numpy.zeros([T.size, 1])

    for T_i in range(T.size):
    for t in range(int(t_end / h + 1)):
    Gamma = numpy.random.normal(size=[n, 1]) * numpy.sqrt(2 * gamma * T[T_i, 0] * h)
    fx1 = p
    fp1 = - x - gamma * p
    x2 = x + fx1 * h
    p2 = p + fp1 * h + Gamma
    fx2 = p2
    fp2 = -x2 - gamma * p2
    x += (fx1 + fx2) / 2 * h
    p += (fp1 + fp2) / 2 * h + Gamma
    H[T_i, 0] = ((x ** 2 + p ** 2) / 2).mean()
    print(H[T_i, 0])

    matplotlib.pyplot.figure()
    matplotlib.pyplot.plot(T, H[:, 0], label="计算值", marker="*", linestyle="")
    matplotlib.pyplot.plot(T, T, label="理论值", linestyle="-")
    matplotlib.pyplot.legend(loc='upper left')
    matplotlib.pyplot.xlabel("T")
    matplotlib.pyplot.ylabel("H")
    matplotlib.pyplot.show()