华中师范大学 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_thickness

return -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 / N

a = 3 * w

return -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 = 0

self.line, = ax.plot([], [], 'k-')

self.x = x

self.ax = ax

# Set up plot parameters

self.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.01

psi = 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 版本的一维薛定谔方程数值解的封装，方便对自定义的势阱计算能级与波函数。

参考文献：

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

pip install SciencePlots

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

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