1. Formalism

为了简化，设一维谐振子的哈密顿量为

\[\hat{H} = - \frac{d^2}{dx^2} + x^2 \]

引入拟设波函数

\[\psi(x, \alpha) = e^{-\alpha x^2}, \]

定义局域能量

\[E_L(x, \alpha) = \frac{ \hat{H} \psi(x, \alpha) }{ \psi(x, \alpha)}, \]

如果 \(\alpha\) 的值恰当，使得 \(\psi(x,\alpha)\) 是 \(\hat{H}\) 的本征态，则 \(\forall x, E_L(x,\alpha) = \epsilon\)。
变分法中的“平均能量”为

\[\langle E(\alpha) \rangle = \int E_L(x, \alpha) \psi(x,\alpha)^2 dx. \]

所以，设置大量的随机行者，用 Markov-Metropolis 方法，使得行者的分布密度 \(\rho(x,\alpha) \propto \psi(x,\alpha)^2\)，然后做 \(E_L(x,\alpha)\)的抽样平均，即得 \(\langle E(\alpha) \rangle\) 的近似值。
另外，可以观察抽样方差 \(\sigma^2(\alpha) = \langle E(\alpha)^2 \rangle - \langle E(\alpha) \rangle^2\) 与 \(\alpha\) 的函数关系。如果 \(\alpha\) 恰当，使得 \(\psi(x,\alpha)\) 是 \(\hat{H}\) 的本征态，则 \(\sigma^2(\alpha) = 0\)，如果 \(\alpha\) 不够理想，就有 \(\sigma^2(\alpha) > 0\)。

2. 代码


# 一维谐振子求基态能量
import numpy as np
import matplotlib.pyplot as plt

# trial wave: psi = e^{ -alpha x^2 }
# H = - d^2/dx^2 + x*x
import scipy.integrate

def psi(x,alpha):
    return np.exp( -alpha * x * x )

#help(scipy.integrate.quad); exit(1)
def plotpsi2(alpha):
    C = scipy.integrate.quad(psi,-np.inf,np.inf,args=(2*alpha,))[0]
    x = np.arange(-3,3,0.1)
    y = [  psi(t,alpha) * psi(t,alpha) / C for t in x ]
    plt.plot(x,y)
#plotpsi2(0.4); exit(1)

# d/dx psi = -2 alpha x e^{- alpha * x * x }
# d^2 /d x^2 psi = -2 alpha + 4 * alpha^2 x^2 e^{-alpha x^2}
# ( - d^2/dx^2 + x*x ) psi / psi =  2 alpha + (- 4 alpha^2 + 1) x^2
def EL( alpha, x):
    return 2*alpha + (1-4*alpha*alpha)*x*x

# if alpha = 0.5, EL(alpha, x) = 1
# if alpha = 0.6, EL(alpha, x) = 1.2 - 0.44 * x*x

def Accept( x1, x2, alpha, psi ):
    if( abs(x2)>3 ): return 0
    p = ( psi(x2,alpha) / psi(x1,alpha) )**2
    #print("p=",p)
    if p>1: return 1
    else: return p

def tryonestep(x,nwalker,h,alpha):
    for j in range(nwalker):
        step = h * ( 2*np.random.randint(0,2) - 1 )
        A = Accept(x[j], x[j] + step, alpha, psi)
        if (np.random.random() < A):
            x[j] += step

def ELsampling( alpha, nwalker, psi, EL, a, b, h, nstep ):
    x = (b-a)*np.random.random(nwalker) + a # 初始位置
    for i in range(nstep):
        tryonestep(x, nwalker, h, alpha)
    #plotpsi2(alpha);
    #plt.hist(x, bins=np.arange(-3.1,3.1,0.2), rwidth=0.8, density="True" ); plt.show()
    #print( "alpha = ",alpha, " <E>_L = ", np.average( [ EL(alpha,t) for t in x ] ) )
    return [ EL(alpha,t) for t in x ]

aveE = []; aveE2 = []　# < E >, < E^2 >
Alpha = np.arange(0.4,0.6,0.02)
for alpha in Alpha:
    ELsample = ELsampling(alpha, 1000, psi, EL, a=-3, b=3, h=0.2, nstep=100 )
    aveE.append( np.average(ELsample) )
    aveE2.append( np.average([t*t for t in ELsample ]) )
sigma2 = aveE2 - np.array([t*t for t in aveE])
plt.plot(Alpha, sigma2)
plt.show()

3. 结果

做出来 \(\sigma^2(\alpha) \sim \alpha\) 的图如下。可见 \(\alpha=0.5\) 是最优参数。

4. 小结

代码中有绘制 \(\rho(x,\alpha)\) 的片段，可以用于观察随机行者在接受概率的指引下收敛于 \(\rho(x,\alpha)\) 的速度。在我的试验中，行者处在 \(-3, -2.8, \cdots, 2.8, 3.0\) 这些格点上，一共大约 30 个位置，大约 100 步随机行走后，行者分布就已经收敛至 \(\rho(x, \alpha)\) 曲线。
\(\langle E(\alpha) \rangle\) 远远没有 \(\sigma^2(\alpha)\) 犀利。上面的代码以及其中的参数跑出来的 \(\langle E(\alpha) \rangle\) 曲线非常磕碜，但是上面展示的 \(\sigma^2(\alpha)\) 曲线就非常规整，且容易反映出问题。如果你开个根号，画 \(\sigma \sim \alpha\)，就会更加 sharp。