用diag直接使用错误_用python学量子力学（1）

华中师范大学 hahakity

在网上看到一个使用 Matlab 教量子力学的文章，很有意思。这里用 python 语言实现一遍，让同学们对量子力学，对偏微分方程的差分近似解法有一个更直观的理解。

学习目标：

理解量子力学的波函数表示与矩阵表示的等价性
学会用向量表示函数，用矩阵表示算符（一阶微分，二阶微分）
学会数值求解任意势阱下定态薛定谔方程的能级与波函数

预备知识：

微分的差分近似
量子力学基础（薛定谔方程）

波函数的向量表示

在用 python 画图时，我们一般先将区间离散化，计算出离散坐标上的函数值，然后画折线图。比如对于函数

, 使用如下代码画图，

# np.linspace 将区间 [-2, 2] 离散化为 100 个坐标点
x = np.linspace(-2, 2, 100)
# 计算 100 个坐标点上的函数值 f
f = np.sin(x) / x
plt.plot(x, f)

将波函数表示为离散坐标点上的实数或复数，写为列向量

变得非常容易理解。

回忆微分的有限差分近似，对于一阶微分，

对于二阶微分，

对区间 [a, b] 所有离散坐标上的 f(x) 微分和二阶微分可以矩阵化，

算符的矩阵表示

当 f(x) 用列向量

表示时，就可以用矩阵来表示微分算子。

因此波动力学与矩阵力学统一。

一阶微分

可以用微分算子矩阵 D 点乘

计算。

二阶微分

可以用微分算子矩阵 D 从左边连续作用两次到

上，也可以使用差分格式直接构造拉普拉斯矩阵，由

计算。

对于随便给定的函数

, 可以看到上述有限差分矩阵作用在离散的

上的结果与解析解非常一致。

解定态薛定谔方程

的任务转化为求哈密顿矩阵 H 的本征值 E 和本征向量

。

注意在量子力学里，动能项里的动量 p 换成了微分算符，进一步用 Laplacian 矩阵表示，势能 U(x) 在区间离散化后，填充在矩阵的对角元上。

如果是多粒子系统，比如考虑多个电子两两之间的库伦相互作用，则两粒子势能

会带来非对角元。

定态薛定谔方程的数值解

这里我用 python 把一维定态薛定谔方程的数值解封装成一个类，后面研究不同势能下的薛定谔方程比较方便。

class Schrodinger:def __init__(self, potential_func, mass = 1, hbar=1,xmin=-5, xmax=5, ninterval=1000):self.x = np.linspace(xmin, xmax, ninterval)    self.U = np.diag(potential_func(self.x), 0)        self.Lap = self.laplacian(ninterval)        self.H = - hbar**2 / (2*mass) * self.Lap + self.U       self.eigE, self.eigV = self.eig_solve()def laplacian(self, N):'''构造二阶微分算子：Laplacian'''dx = self.x[1] - self.x[0]return (-2 * np.diag(np.ones((N), np.float32), 0)+ np.diag(np.ones((N-1), np.float32), 1)+ np.diag(np.ones((N-1), np.float32), -1))/(dx**2)def eig_solve(self):'''解哈密顿矩阵的本征值，本征向量；并对本征向量排序'''w, v = np.linalg.eig(self.H)  idx_sorted = np.argsort(w)     return w[idx_sorted], v[:, idx_sorted]def wave_func(self, n=0):return self.eigV[:, n]def eigen_value(self, n=0):return self.eigE[n]def check_eigen(self, n=7):'''check wheter H|psi> = E |psi> '''with plt.style.context(['science', 'ieee']):HPsi = np.dot(self.H, self.eigV[:, n])EPsi = self.eigE[n] * self.eigV[:, n]plt.plot(self.x, HPsi, label=r'$H|psi_{%s} rangle$'%n)plt.plot(self.x, EPsi, '-.', label=r'$E |psi_{%s} rangle$'%n)plt.legend(loc='upper center')plt.xlabel(r'$x$')plt.ylim(EPsi.min(), EPsi.max() * 1.6)def plot_density(self, n=7):with plt.style.context(['science', 'ieee']):rho = self.eigV[:, n] * self.eigV[:, n]plt.plot(self.x, rho)plt.title(r'$E_{%s}=%.2f$'%(n, self.eigE[n]))plt.ylabel(r'$rho_{%s}(x)=psi_{%s}^*(x)psi_{%s}(x)$'%(n, n, n))plt.xlabel(r'$x$')def plot_potential(self):with plt.style.context(['science', 'ieee']):plt.plot(self.x, np.diag(self.U))plt.ylabel(r'potential')plt.xlabel(r'$x$')

谐振子势

# 定义谐振子势
def harmonic_potential(x, k=100):return 0.5 * k * x**2# 创建谐振子势下的薛定谔方程
schro_harmonic = Schrodinger(harmonic_potential)

从上面例子可以看到，封装的比较完整，对任意 1 维势调用很简单。先来可视化谐振子势能，

schro_harmonic.plot_potential()

schro_harmonic.check_eigen(n=1)

再用上面这条命令检查一下薛定谔方程的解是否准确，具体来说就是本征方程

是否满足。

还可以看看粒子在谐振子势阱中的分布概率密度，

# 这里随便选了一个能级 n = 9
schro_harmonic.plot_density(n=9)

如果亲自尝试一下，你会发现在上面这个谐振子势下，数值解的能级与解析解非常接近

对比解析解，

其中

。

解析解中，n=0 时，

。 n = 1, 2, 3... 时，

是

的 3 倍，5倍，7倍... 。

Woods Saxon 势能

在核物理领域，原子核中一大团核子所产生的势能接近于 Woods Saxon 函数形式。这里看看Woods Saxon 势阱中一个核子的能级分布。势阱函数形式为，

def woods_saxon_potential(x, R0=6.2, surface_thickness=0.5):sigma = surface_thicknessreturn  -1 / (1 + np.exp((np.abs(x) - R0)/sigma))

用这个势阱构造薛定谔方程，

ws_schro = Schrodinger(woods_saxon_potential)

先画一下势阱的样子，

ws_schro.plot_potential()

再画一下 Woods Saxon 势阱中核子的波函数和能级，

对于基态 n=0

ws_schro.plot_density(n=0)

第 n=9 激发态

与谐振子势阱的结果有相当大差别。

双势阱

def double_well(x, xmax=5, N=100):w = xmax / Na = 3 * wreturn -100 * (np.heaviside(x + w - a, 0.5) - np.heaviside(x - w - a, 0.5)+np.heaviside(x + w + a, 0.5) - np.heaviside(x - w + a, 0.5))dw = lambda x: double_well(x, xmax=5, N=1000)
dw_shro = Schrodinger(double_well)

双势阱中前几个能级下粒子的概率密度分布，

dw_shro.plot_density(n=0)
dw_shro.plot_density(n=1)
dw_shro.plot_density(n=2)
dw_shro.plot_density(n=4)

双势阱中粒子概率密度分布的随时间演化

下面这个问题考虑一个粒子被捕获在上例所示的有限深方势阱中，初态为基态与第一激发态的叠加态，观察粒子的概率密度分布随时间的演化。初态为，

这里画图看一下基态

，第一激发态

和两者叠加态

的波函数，

psi0 = dw_shro.wave_func(n=0)
psi1 = dw_shro.wave_func(n=1)
psi = 1 / np.sqrt(2) * (psi0 + psi1)with plt.style.context(['science', 'ieee']):plt.plot(dw_shro.x, psi0, 'r--', label=r'$|Psi_{E_0} rangle$')plt.plot(dw_shro.x, psi1, 'b:', label=r'$|Psi_{E_1} rangle$')plt.plot(dw_shro.x, psi, 'k-', label=r'$|Psi(t=0)rangle = (|Psi_{E_0}rangle + |Psi_{E_1}rangle) / sqrt{2}$')plt.legend(loc='best')plt.xlabel(r'$x$')plt.ylim(-0.3, 0.3)plt.xlim(-2, 2)

如下图所示，初始时刻基态与第一激发态在右边势阱处相消，导致叠加态的波函数在左边势阱处有峰值结构。后面演示此峰如何随时间在两个势阱间振荡。

叠加态波函数的时间演化直接用时间演化算符，

def psit(t, hbar=1):'''基态与第一激发态的叠加态波函数，随时演化'''psi0 = dw_shro.wave_func(n=0)psi1 = dw_shro.wave_func(n=1)E0 = dw_shro.eigen_value(0)E1 = dw_shro.eigen_value(1)return 1/np.sqrt(2) * (psi0 * np.exp(-1j * E0 * t/hbar)+  psi1 * np.exp(-1j * E1 * t/hbar))

注意我们用 Dirac

列向量表示所有离散空间点上的波函数值。用

表示给定

点上的波函数值。波函数的平方表示概率密度，对于给定的离散时空点

当

函数项的值从 1 变到负 1 ，

点的概率密度从

变到

下面是概率密度在两个势阱间震荡间振荡的动画，

知乎视频www.zhihu.com

# 动画代码
%matplotlib notebook
from matplotlib.animation import FuncAnimation
class UpdateDist:def __init__(self, ax, x):self.success = 0self.line, = ax.plot([], [], 'k-')self.x = xself.ax = ax# Set up plot parametersself.ax.set_xlim(-0.6, 0.6)self.ax.set_ylim(-0.02, 0.1)self.ax.grid(True)def __call__(self, i):time = i * 0.01psi = psit(t = time)density = np.real(np.conjugate(psi) * psi)     self.line.set_data(self.x, density)return self.line,# 画缩放了的双势阱
potential = double_well(dw_shro.x) * 1.0E-4
fig, ax = plt.subplots()
ax.plot(dw_shro.x, potential, ':')
ax.set_xlabel(r'$x$')
ud = UpdateDist(ax, x=dw_shro.x)
anim = FuncAnimation(fig, ud, frames=1000, interval=100, blit=True)
#anim.save('../htmls/images/double_well_evolution.mp4')
plt.show()

总结：

使用微分的有限差分近似可以将波函数表示为向量，将微分算子化为矩阵，将定态薛定谔方程的求解化为哈密顿矩阵的本征值，本征向量求解问题。实现了 python 版本的一维薛定谔方程数值解的封装，方便对自定义的势阱计算能级与波函数。

参考文献：

【1】https://arxiv.org/pdf/0704.1622.pdf

【2】Teaching Quantum Mechanics with MATLAB

注：画图用到了 matplotlib 库，要得到本文画图风格，需要安装 SciencePlots 库。

pip install SciencePlots

如果出错，开启 no-latex 选项。

with plt.style.context(["science", "ieee", "no-latex"]):