Sympy中量子谐振子哈密顿量矩阵元素结果不符合预期

一、非对角元异常（i+2=j时不为0）的原因及修复 #

你的代码中哈密顿量的势能项写错了：标准量子谐振子的势能是$\frac{1}{2}m\omega^2x2$，但你写成了1/2*m**2*omega**2*x**2（多了一个m的平方）。这个错误导致哈密顿量偏离标准形式，进而使得j=i+2时，x²的矩阵元贡献出非零值。

修复方法 #

修改哈密顿量定义中的势能项，将m**2改为m：

def H(psi):
    return -hbar**2/(2*m) * sp.diff(psi, x, x) + 1/2*m*omega**2*x**2*psi

修正后，标准谐振子哈密顿量的非对角元会全部为0，对角元为$(n+\frac{1}{2})\hbar\omega$，符合预期。

二、hbar**2/hbar未化简的问题解决 #

你使用的sympy.physics.units.hbar是带物理单位的量，SymPy对带单位的幂次表达式化简需要特殊处理，有两种解决方案：

方案1：使用无单位的hbar符号 #

直接将hbar定义为普通符号，而非带单位的量，SymPy会自动完成幂次化简：

hbar = sp.Symbol('hbar')

方案2：保留单位时强制化简 #

如果需要保留物理单位，可以在计算结果后用sp.powsimp()强制处理幂次项：

# 在matrix_elem函数中修改返回值
return sp.powsimp(sp.simplify(result))

例如对hbar**2/hbar执行sp.powsimp()会直接得到hbar。

修正后的完整代码 #

import sympy as sp
from sympy.physics.qho_1d import psi_n
import numpy as np

x = sp.Symbol('x', real=True)
m = 1
omega = 1
# 选择方案1：无单位hbar，自动化简
hbar = sp.Symbol('hbar')
# 若选方案2，保留单位：
# hbar = sp.physics.units.hbar

def H(psi):
    # 修正后的哈密顿量
    return -hbar**2/(2*m) * sp.diff(psi, x, x) + 1/2*m*omega**2*x**2*psi

def matrix_elem(H, psi1, psi2):
    integrand = psi1.conjugate() * H(psi2)
    integral_result = sp.integrate(sp.simplify(integrand), (x, -sp.oo, sp.oo))
    # 方案2需添加powsimp
    return sp.powsimp(sp.simplify(integral_result))

max_n = 5
H_matrix = sp.Matrix(
    [[matrix_elem(H, psi_n(i, x, m, omega), psi_n(j, x, m, omega)) for j in range(max_n)] for i in range(max_n)]
)

print(sp.simplify(H_matrix))
# 检查非对角元是否为0
for i in range(max_n):
    for j in range(max_n):
        if i != j and H_matrix[i,j] != 0:
            print(f"{i} {j}: {H_matrix[i,j]=}")

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