Scipy求解一维氢原子哈密顿量本征态异常原因排查

问题背景 #

我在求解一维氢原子哈密顿量本征态时，按定义空间网格、构造哈密顿量、使用scipy.linalg.eigh计算本征值与本征向量的步骤实现，但结果出现异常：本征向量符合量子理论预期，但半数本征值存在错误——部分理论上不该存在的本征值出现，部分应存在的本征值错误关联了对应的本征向量（单位为原子单位）。

物理分析 #

• 能量（本征值）与波函数（本征向量）需符合量子理论；
• 我原本认为一维氢原子束缚态能量满足 ( E_n = -\frac{0.5}{n^2} )，据此判断计算结果中第1、3、5个本征值不应存在；
• 橙色曲线的波函数形态符合 ( E_n = -\frac{0.5}{2^2} )（即n=2）的理论波函数，但对应能量值错误；
• 首个本征值 ( E_1 ) 的取值受epsilon参数与网格点数影响极大。

所用代码 #

import numpy as np
from scipy.sparse import diags
from scipy.linalg import eigh

import matplotlib.pyplot as plt

# Parameters
L = 50.0         # domain size
N = 1000         # number of grid points
x = np.linspace(-L/2, L/2, N)
dx = x[1] - x[0]

# 1D Hydrogen-like potential (softened to avoid singularity at x=0)
Z = 1.0
epsilon = 1e-2
V = -Z / np.sqrt(x**2 + epsilon**2)

# Kinetic energy operator (finite difference, central)
main_diag = np.full(N, -2.0)
off_diag = np.full(N-1, 1.0)
laplacian = diags([off_diag, main_diag, off_diag], [-1, 0, 1]) / dx**2

# Hamiltonian
H = -0.5 * laplacian.toarray() + np.diag(V)

# Solve eigenvalue problem
num_states = 10
eigvals, eigvecs = eigh(H)
eigvals = eigvals[:num_states]
eigvecs = eigvecs[:, :num_states]

# Normalize eigenfunctions
eigvecs /= np.sqrt(np.sum(np.abs(eigvecs)**2, axis=0) * dx)

# Plot
plt.figure(figsize=(8,6))
for i in range(num_states):
    plt.plot(x, eigvecs[:, i] + eigvals[i], label=f'n={i+1}, E={eigvals[i]:.3f}')
plt.xlabel('x')
plt.ylabel('Wavefunction + Energy')
plt.title('1D Hydrogen Atom Eigenstates')
plt.legend()
plt.grid()

plt.figure()
plt.plot(x, eigvecs[:, 0])
plt.plot(x, np.exp(-np.abs(x)), label='Analytical 1s state', linestyle='--', color='red')
plt.show()

异常原因分析 #

1. 错误套用三维氢原子能级公式 #

你混淆了一维与三维氢原子的能级特性：

• 三维氢原子的能级公式 ( E_n = -\frac{0.5}{n^2} ) 并不适用于一维系统；
• 一维氢原子的束缚态能量为 ( E_k = -\frac{0.5}{(k + 0.5)^2} )（原子单位），其中 ( k = 0,1,2,... )；
• 每个能级 ( E_k ) 对应两个宇称简并态（偶宇称和奇宇称），因此计算结果中会出现成对的相近本征值，这是一维系统的固有特性，并非错误。你看到的“不该存在”的本征值，其实是同一能级的简并态。

2. 奇点软化参数epsilon的干扰 #

为避免x=0处的奇点，你引入了epsilon参数，但：

• 当epsilon取值过大（如1e-2）时，势阱的深度和形状会被显著修改，导致基态能量偏离理论值；
• 网格点数不足时，空间分辨率无法准确捕捉波函数在奇点附近的急剧变化，进一步放大了能量计算的误差。

3. 边界条件与网格范围的偏差 #

你选择的网格范围 ( x \in [-25,25] ) 可能不足以让波函数充分衰减：

• 有限差分法隐含了Dirichlet边界条件（波函数在网格边界处为0），若波函数在边界处仍有可观振幅，会导致能级计算出现偏差；
• 较小的网格范围会人为截断波函数的尾部，影响能量本征值的准确性。

修正建议 #

1. 采用正确的一维能级公式验证：用 ( E_k = -\frac{0.5}{(k+0.5)^2} ) 核对计算结果，而非三维氢原子的公式；
2. 优化奇点处理：将epsilon减小至1e-4或更小，同时增加网格点数N（如2000），平衡奇点软化与数值精度；
3. 扩大网格范围：将L增大至100，确保波函数在边界处的振幅衰减到可忽略的程度；
4. 使用稀疏矩阵求解：避免将稀疏拉普拉斯矩阵转为稠密矩阵，改用scipy.sparse.linalg.eigsh直接求解稀疏哈密顿量，提升计算效率与精度。

内容的提问来源于stack exchange，提问作者Maxence BARRÉ