Wigner方程 particle random walk随机粒子方法 代码

Wigner方程描述了量子力学中粒子的运动，是一个偏微分方程。而随机粒子方法是一种数值模拟方法，可以用来模拟Wigner方程的解，特别适用于高维量子系统的数值计算。下面给出使用随机粒子方法模拟Wigner方程的代码示例。

首先，我们需要导入必要的库：

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

然后，我们定义一些常数和初始条件：

hbar = 1.0
m = 1.0
omega = 1.0
sigma = 1.0
x0 = 0.0
p0 = 1.0
Nt = 1000
dt = 0.1
Np = 1000

其中，hbar、m、omega、sigma分别表示普朗克常数、质量、频率和初始粒子分布的标准差，x0和p0为初始位置和动量，Nt为时间步数，dt为时间步长，Np为随机粒子数。

接下来，我们定义一个函数，用于计算Wigner函数的演化：

def wigner_evolution(x, p, t):
    def V(x):
        return 0.5 * m * omega**2 * x**2
    V_x = V(x)
    V_p = 0.0
    dx = np.sqrt(2.0 * dt / hbar) * np.random.normal(size=Np)
    dp = np.sqrt(2.0 * dt / hbar) * np.random.normal(size=Np)
    x_new = x + p / m * dt + dx
    p_new = p - V_x / hbar * dt - p / m * dt + dp
    W_new = np.zeros((Np, Np))
    for i in range(Np):
        for j in range(Np):
            W_new[i, j] = np.exp(
                -0.5 * ((x_new[i] - x_new[j])**2 / sigma**2
                        + (p_new[i] - p_new[j])**2 * sigma**2)) / (2.0 * np.pi * hbar)
    return W_new

其中，V(x)为势函数，dx和dp分别为随机粒子的位置和动量扰动，x_new和p_new为演化后的位置和动量，W_new为演化后的Wigner函数。

最后，我们可以使用Euler方法来数值求解Wigner方程，并画出Wigner函数的演化：

x = np.random.normal(loc=x0, scale=sigma, size=Np)
p = np.random.normal(loc=p0, scale=sigma, size=Np)
W = np.zeros((Np, Np))
for i in range(Np):
    for j in range(Np):
        W[i, j] = np.exp(
            -0.5 * ((x[i] - x[j])**2 / sigma**2
                    + (p[i] - p[j])**2 * sigma**2)) / (2.0 * np.pi * hbar)
for it in range(Nt):
    W = wigner_evolution(x, p, it*dt)
    x = x + p / m * dt
    p = p - integrate.quad(lambda x: V(x), x - 0.5 * sigma, x + 0.5 * sigma)[0] / hbar * dt - p / m * dt
    W = np.mean(W, axis=0)
    plt.plot(x, W)
plt.show()

运行以上代码即可得到Wigner函数的演化图像。注意，由于随机粒子方法本身的随机性，每次运行结果可能会略有不同。