华中师范大学 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
# 动画代码
%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"]):
更多推荐
python求解薛定谔方程_用python学量子力学(1)
发布评论