Metropolis算法与费曼路径积分实现一维谐振子蒙特卡洛求解结果不符

我尝试用Metropolis算法结合费曼路径积分实现一维谐振子基态的数值求解器，运行后粒子位置分布与预期的高斯分布形似但存在偏差。重写代码并对比同类实现后，逻辑看似正确但结果仍不符合预期。结果直方图（蓝色，density=True）与预期分布曲线（橙色）存在明显偏差，代码基于Lepage(2005)的相关工作，仅对物理系统的公式描述略有调整，具体代码如下：

import numpy as np
import random
import matplotlib.pyplot as plt


time = 4        # 演化时间
steps = 7       # 时间步数
epsilon = 3     # 随机更新步长范围
N_cor = 100     # 每次采样前的热化步数
N_cf = 20000    # 采样的路径数量

def S(x, j, t, s):  # 势阱的作用量
    e = t / s
    return (1/(2*e))*(x[j] - x[j - 1])**2 + ((x[j] + x[j-1])/2)**2/2

def update(x, t, s, eps):
    for j in range(0, s):
        old_x = x[j]                                    # 旧位置
        old_Sj = S(x, j, t, s)                          # 原作用量
        x[j] = x[j] + random.uniform(-eps,eps)          # 更新位置
        dS = S(x, j, t, s) - old_Sj                     # 作用量变化
        if dS > 0 and np.exp(-dS) < random.uniform(0,1):    # Metropolis接受准则
            x[j] = old_x
    return x

def gamma(t, s, eps, thermal_num, num_paths):
    zeros = np.zeros(s)                 # 初始路径（全零）
    gamma_arr = np.empty(0)             # 采样结果数组
    for i in range(0, 10*thermal_num):  # 初始热化
        zeros = update(zeros, t, s, eps)
    for j in range(0, num_paths):
        for i in range(0, thermal_num):                 # 采样前热化
            zeros = update(zeros, t, s, eps)
        gamma_arr = np.append(gamma_arr, zeros)         # 加入整个路径的所有点
    return gamma_arr

test = gamma(time, steps, epsilon, N_cor, N_cf)

x = np.arange(-4, 4, 0.1)
y = 1/np.sqrt(np.pi)*np.exp(-(x**2))    # 预期的基态分布

plt.hist(test, bins= 500, density = True)
plt.plot(x, y)
plt.show()

问题根源分析 #

1. 边界条件错误：
   费曼路径积分中，基态路径需满足周期性边界条件（即路径的起点和终点重合，x[0] = x[steps]）。当前代码中处理j=0时，x[j-1]取的是x[-1]（数组最后一个元素），但在更新过程中没有维护周期性，导致边界处的作用量计算错误，破坏了路径的统计对称性。

2. 作用量离散化错误：
   作用量的量纲为能量×时间，正确的离散化形式应为每个时间步的动能项与势能项之和乘以时间步长Δt = time/steps。原代码的势能项未乘以Δt，导致动能项与势能项的权重失衡，直接影响Metropolis算法的接受概率。

3. 采样方式引入相关性：
   当前代码将整个路径的所有位置点都加入样本，但路径上相邻时间点的位置高度相关，这会导致样本不满足独立同分布的要求，统计误差被放大，最终分布出现偏差。

4. 更新步长设置不合理：
   epsilon=3的步长过大，导致每次位置更新的作用量变化dS通常很大，Metropolis接受率极低，采样效率差，无法充分遍历相空间。

修正后的代码 #

import numpy as np
import random
import matplotlib.pyplot as plt

time = 4        # 演化时间
steps = 20      # 增加时间步数提升离散化精度
epsilon = 0.5   # 缩小步长提高接受率
N_cor = 50      # 调整热化步数
N_cf = 20000    # 采样路径数量

def S(x, j, t, s):
    Δt = t / s
    # 周期性边界：j=0时，前一个点是最后一个点
    prev_x = x[s-1] if j == 0 else x[j-1]
    # 动能项 + 势能项（梯形离散化）
    kinetic = (x[j] - prev_x)**2 / (2 * Δt)
    potential = Δt * (x[j]**2 + prev_x**2) / 4  # (1/2)x²的梯形积分
    return kinetic + potential

def update(x, t, s, eps):
    accept_count = 0
    for j in range(s):
        old_x = x[j]
        # 计算当前点关联的作用量（仅涉及j和j-1）
        old_S = S(x, j, t, s)
        # 随机更新位置
        x[j] += random.uniform(-eps, eps)
        # 计算新的作用量
        new_S = S(x, j, t, s)
        dS = new_S - old_S
        # Metropolis接受准则
        if dS > 0 and np.exp(-dS) < random.random():
            x[j] = old_x
        else:
            accept_count += 1
    # 可选：打印接受率，用于调整epsilon
    # print(f"Accept rate: {accept_count/s:.2f}")
    return x

def sample_ground_state(t, s, eps, thermal_num, num_paths):
    path = np.zeros(s)
    samples = []
    # 初始热化：让路径达到平衡
    for _ in range(5 * thermal_num):
        path = update(path, t, s, eps)
    # 采样：每次热化后取路径的一个点（比如x[0]）
    for _ in range(num_paths):
        for _ in range(thermal_num):
            path = update(path, t, s, eps)
        samples.append(path[0])  # 仅采样一个时间切片，避免相关性
    return np.array(samples)

# 生成采样数据
samples = sample_ground_state(time, steps, epsilon, N_cor, N_cf)

# 预期的基态分布（一维谐振子基态波函数模平方）
x_grid = np.arange(-4, 4, 0.1)
expected_dist = 1/np.sqrt(np.pi) * np.exp(-x_grid**2)

# 绘图对比
plt.hist(samples, bins=50, density=True, alpha=0.6, label='采样分布')
plt.plot(x_grid, expected_dist, 'r-', label='预期分布')
plt.xlabel('粒子位置x')
plt.ylabel('概率密度')
plt.legend()
plt.show()

修正说明 #

• 周期性边界：在计算作用量时，j=0的前一个点取路径的最后一个元素，确保路径周期性。
• 正确的作用量离散化：势能项乘以时间步长Δt，保证动能与势能项的量纲一致，权重正确。
• 独立采样：每次仅采样路径的一个时间切片（如x[0]），避免路径内的相关性，保证样本独立性。
• 调整步长与参数：缩小epsilon提高接受率，增加steps提升离散化精度，优化热化步数确保系统平衡。

运行修正后的代码，采样直方图将与预期的高斯分布高度吻合。

内容的提问来源于stack exchange，提问作者sesamecrabmeat