Aubry-André哈密顿量波矢时间演化数值复现问题求助

核心错误：索引混淆导致方差计算错误 #

你代码中方差计算的核心问题是用时间循环的索引n替代了格点位置索引。原公式$var = \sqrt{\sum(n - n_0)^2 |\psi(t)|^2}$中的n是晶格的格点位置序号，而不是时间步的循环下标，这直接导致计算结果完全偏离预期。

修正后方差计算代码 #

首先提前定义与波函数分量对应的格点位置数组（放在n0定义之后）：

# 定义格点位置数组，对应psi_0每个分量的晶格位置
# 若你的晶格是1到N-1（共N-1个格点），则用：
grid = np.arange(1, N)
# 若直接用数组索引作为格点位置，则用：
# grid = np.arange(N-1)

然后修改方差计算的循环逻辑，替换错误的索引：

variance = []
for psi_evo in psi_t:
    prob_density = np.abs(psi_evo) ** 2
    # 使用格点位置数组grid计算平方位移的加权和
    mean_sq_displacement = np.sum((grid - n0) ** 2 * prob_density)
    variance.append(np.sqrt(mean_sq_displacement))

额外优化与注意点 #

1. 初始态索引匹配：当N=100时，n0=(N+1)//2=50，但psi_0是长度为N-1=99的数组（索引0到98），若你的晶格位置是1到99，那么中间位置50对应的数组索引应为49，需确认初始态psi_0[n0]=1是否对应正确的晶格中间位置。
2. 时间演化效率优化：当前代码中每次循环都构造对角矩阵U，可以改为直接对本征态分量乘相位因子，避免重复构造矩阵，提升效率：

# 替换原时间演化循环
psi_t = []
for t in tlim:
    # 直接计算每个本征态分量的相位演化
    phi_t = phi_0 * np.exp(-1j * energias * t)
    psi_evo = np.dot(S, phi_t)
    psi_t.append(psi_evo)

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