有限差分求解定态

正文里能精确求解的势场屈指可数。但只要肯用电脑，任何一维势场的束缚态 都可以在一秒钟内解出来——方法就是把微分方程变成矩阵本征值问题， 交给成熟的对角化例程。4.4 节与 4.5 节的自定义势求解器用的正是这套算法， 这里给出完整的推导、代码与检验。

思路：函数变向量，算符变矩阵

定态方程（自然单位 $\hbar=m=1$，代码里统一用这套单位）：

$$-\frac{1}{2}\psi''(x)+V(x)\psi(x)=E\psi(x)\tag{C.1.1}$$

离散化：把连续区间切成 $N$ 个格点 $x_j$，间距 $h$。 波函数 $\psi(x)$ 变成一个 $N$ 维向量 $(\psi_1,\dots,\psi_N)$—— 这正是第 3 章「函数就是无穷维向量」 的字面实现，只不过现在维数有限。

剩下的问题：二阶导数 $\psi''$ 怎么用格点值表示？

三点差分格式与它的误差阶基础~5 min展开

把 $\psi(x_j\pm h)$ 在 $x_j$ 处 Taylor 展开：

$$\psi_{j\pm1}=\psi_j\pm h\psi_j'+\frac{h^2}{2}\psi_j''\pm\frac{h^3}{6}\psi_j'''+\frac{h^4}{24}\psi_j^{(4)}+\cdots\tag{C.1.2}$$

两式相加：奇数阶项（$\psi'$、$\psi'''$）成对相消，

$$\psi_{j+1}+\psi_{j-1}=2\psi_j+h^2\psi_j''+\frac{h^4}{12}\psi_j^{(4)}+\cdots\tag{C.1.3}$$

整理得三点差分格式：

$$\psi_j''=\frac{\psi_{j+1}-2\psi_j+\psi_{j-1}}{h^2}-\frac{h^2}{12}\psi_j^{(4)}+\cdots\tag{C.1.4}$$

误差首项 $\propto h^2$：格距减半，误差变成四分之一。这叫二阶精度， 待会儿我们用数值实验直接验证这个「4 倍」。

代入定态方程，第 $j$ 行变成只连接相邻格点的代数方程：

$$-\frac{1}{2h^2}\psi_{j-1}+\left(\frac{1}{h^2}+V_j\right)\psi_j-\frac{1}{2h^2}\psi_{j+1}=E\psi_j\tag{C.1.5}$$

写成矩阵形式，$H$ 是三对角矩阵：主对角元 $1/h^2+V_j$， 两条次对角元全是 $-1/2h^2$。边界条件：取足够大的盒子并要求 $\psi_0=\psi_{N+1}=0$（Dirichlet 边界，相当于套一个无限深势阱）—— 只要盒子比波函数的衰减范围大得多，这个人为的墙就没有影响。

完整代码

解谐振子 $V=x^2/2$，与精确解 $E_n=n+\tfrac12$ 对比：

import numpy as np
from scipy.linalg import eigh_tridiagonal

# ---- 参数（自然单位 hbar = m = 1）----
L = 16.0          # 求解区间 [-L/2, L/2]
N = 2000          # 内部格点数
x = np.linspace(-L/2, L/2, N + 2)[1:-1]   # 去掉两端边界点
h = x[1] - x[0]

# ---- 势能：谐振子 V = x^2 / 2（换势场只改这一行）----
V = 0.5 * x**2

# ---- 哈密顿矩阵：三对角 ----
diag = 1.0 / h**2 + V                     # 主对角元
off  = -0.5 / h**2 * np.ones(N - 1)       # 次对角元

# ---- 对角化（只要最低几个本征值就够了）----
E, psi = eigh_tridiagonal(diag, off, select='i', select_range=(0, 5))

# ---- 归一化：eigh 返回的向量满足 sum(psi^2)=1，换成积分归一 ----
psi /= np.sqrt(h)

print("  n   数值 E_n     精确 E_n     误差")
for n in range(6):
    exact = n + 0.5
    print(f"  {n}   {E[n]:.8f}   {exact:.1f}          {abs(E[n]-exact):.2e}")

运行输出（实测）：

  n   数值 E_n     精确 E_n     误差
  0   0.49999800   0.5          2.00e-06
  1   1.49999001   1.5          9.99e-06
  2   2.49997403   2.5          2.60e-05
  3   3.49995005   3.5          5.00e-05
  4   4.49991808   4.5          8.19e-05
  5   5.49987812   5.5          1.22e-04

六位有效数字的精度，耗时不到一秒。逐行说明关键点：

• np.linspace(...)[1:-1]：生成含两端的 $N+2$ 个点后扔掉两端。 边界上 $\psi=0$ 是已知的，不该作为未知数；矩阵里也就自然不含它们—— 这就是 Dirichlet 边界的全部实现。
• eigh_tridiagonal：专吃三对角矩阵的对角化例程，复杂度远低于对满矩阵硬来 （select_range=(0, 5) 表示只算指标 0 到 5 的最低六个本征对，更快）。
• psi /= np.sqrt(h)：数值例程按向量内积 $\sum_j\psi_j^2=1$ 归一， 物理要的是 $\int\psi^2\dd x\approx h\sum_j\psi_j^2=1$，差一个 $\sqrt h$。
• 注意能级越高误差越大：高激发态波函数振荡快，同样的格距分辨得越吃力。

收敛性检验：眼见为实的 h² 定律

任何数值结果都应该做这一步：网格加密一倍，看误差缩小几倍。

import numpy as np
from scipy.linalg import eigh_tridiagonal

def ground_energy(N, L=16.0):
    x = np.linspace(-L/2, L/2, N + 2)[1:-1]
    h = x[1] - x[0]
    diag = 1.0 / h**2 + 0.5 * x**2
    off = -0.5 / h**2 * np.ones(N - 1)
    E, _ = eigh_tridiagonal(diag, off, select='i', select_range=(0, 0))
    return E[0], h

print("   N        h        |E0 - 0.5|    误差比")
prev = None
for N in [125, 250, 500, 1000, 2000]:
    E0, h = ground_energy(N)
    err = abs(E0 - 0.5)
    ratio = f"{prev/err:.2f}" if prev else "  --"
    print(f"{N:5d}   {h:.4f}   {err:.3e}   {ratio}")
    prev = err

实测输出：

   N        h        |E0 - 0.5|    误差比
  125   0.1270   5.044e-04     --
  250   0.0637   1.270e-04   3.97
  500   0.0319   3.187e-05   3.98
 1000   0.0160   7.984e-06   3.99
 2000   0.0080   1.998e-06   4.00

误差比稳定收敛到 4——正是推导里预言的 $O(h^2)$。 如果你的程序跑出来的比值不是 4（比如是 2，或者乱跳）， 说明有别的误差源在捣乱，最常见的就是下面这几个坑。

换一个势场玩玩

代码里 V = 0.5 * x**2 一行换掉即可。几个值得一试的：

势场	代码	看什么
四次振子 $x^4$	V = x**4	无精确解！能级间距随 $n$ 变宽（对比谐振子等间距）
双阱	V = (x**2 - 4)**2 / 8	最低两条能级几乎简并——隧穿劈裂，联系 11.4 瞬子
有限深势阱	V = np.where(abs(x) < 2, -1.0, 0.0)	束缚态个数有限，对比 2.8 节
线性势 $\lvert x\rvert$	V = np.abs(x)	夸克禁闭的玩具模型；能级由 Airy 函数零点给出

∎本节关键公式

• 三点差分

  $$\psi_j''=\frac{\psi_{j+1}-2\psi_j+\psi_{j-1}}{h^2}+O(h^2)$$

  Taylor 展开相加，奇数阶自动相消

• 离散哈密顿量

  $$H_{jj}=\frac{1}{h^2}+V_j,\qquad H_{j,j\pm1}=-\frac{1}{2h^2}$$

  三对角矩阵；Dirichlet 边界 = 直接去掉边界未知数

• 收敛判据

  $$\frac{\text{err}(h)}{\text{err}(h/2)}\to 4$$

  二阶方法的指纹；比值不对说明另有误差源

?自测共 3 题

1. 1.

   把格点间距 h 缩小为原来的 1/3，三点差分的离散化误差大约变为？

   • 1/9
   • 1/3
   • 1/27
   • 不变

   提交
2. 2.

   数值解谐振子时发现：网格已经很密，但基态能量仍偏离 0.5 且不再随 N 改善。最可能的原因是？

   • 盒子长度 L 不够大，波函数被边界的人为墙压住了
   • scipy 的对角化例程精度不足
   • 谐振子本来就没有精确解
   • 忘了把波函数归一化

   提交
3. 3.

   对角化 N = 2000 的哈密顿矩阵会得到 2000 个本征值，其中可信的是？

   • 低端那些波长远大于 h、且在边界前已衰减为零的束缚态
   • 全部——它们都是矩阵的严格本征值
   • 只有基态
   • 高端的，因为高能态受势场细节影响小

   提交