華中師范大學(xué) hahakity

在網(wǎng)上看到一個(gè)使用 Matlab 教量子力學(xué)的文章，很有意思。這里用 python 語(yǔ)言實(shí)現(xiàn)一遍， 讓同學(xué)們對(duì)量子力學(xué)，對(duì)偏微分方程的差分近似解法有一個(gè)更直觀的理解。

學(xué)習(xí)目標(biāo)：

理解量子力學(xué)的波函數(shù)表示與矩陣表示的等價(jià)性

學(xué)會(huì)用向量表示函數(shù)，用矩陣表示算符（一階微分，二階微分）

學(xué)會(huì)數(shù)值求解任意勢(shì)阱下定態(tài)薛定諤方程的能級(jí)與波函數(shù)

預(yù)備知識(shí)：

微分的差分近似

量子力學(xué)基礎(chǔ)（薛定諤方程）

波函數(shù)的向量表示

在用 python 畫(huà)圖時(shí)，我們一般先將區(qū)間離散化，計(jì)算出離散坐標(biāo)上的函數(shù)值，然后畫(huà)折線圖。比如對(duì)于函數(shù)

, 使用如下代碼畫(huà)圖，# np.linspace 將區(qū)間 [-2, 2] 離散化為 100 個(gè)坐標(biāo)點(diǎn) x = np.linspace(-2, 2, 100) # 計(jì)算 100 個(gè)坐標(biāo)點(diǎn)上的函數(shù)值 f f = np.sin(x) / x plt.plot(x, f)

將波函數(shù)表示為離散坐標(biāo)點(diǎn)上的實(shí)數(shù)或復(fù)數(shù)，寫(xiě)為列向量

變得非常容易理解。
回憶微分的有限差分近似，對(duì)于一階微分，

對(duì)于二階微分，

對(duì)區(qū)間 [a, b] 所有離散坐標(biāo)上的 f(x) 微分和二階微分可以矩陣化，

算符的矩陣表示

當(dāng) f(x) 用列向量

表示時(shí)，就可以用矩陣來(lái)表示微分算子。

因此波動(dòng)力學(xué)與矩陣力學(xué)統(tǒng)一。

一階微分

可以用微分算子矩陣 D 點(diǎn)乘 計(jì)算。

二階微分

可以用微分算子矩陣 D 從左邊連續(xù)作用兩次到 上，也可以使用差分格式直接構(gòu)造拉普拉斯矩陣，由 計(jì)算。

對(duì)于隨便給定的函數(shù)

, 可以看到上述有限差分矩陣作用在離散的 上的結(jié)果與解析解非常一致。

解定態(tài)薛定諤方程

的任務(wù)轉(zhuǎn)化為求哈密頓矩陣 H 的本征值 E 和本征向量 。

注意在量子力學(xué)里，動(dòng)能項(xiàng)里的動(dòng)量 p 換成了微分算符，進(jìn)一步用 Laplacian 矩陣表示，勢(shì)能 U(x) 在區(qū)間離散化后，填充在矩陣的對(duì)角元上。

如果是多粒子系統(tǒng)，比如考慮多個(gè)電子兩兩之間的庫(kù)倫相互作用，則兩粒子勢(shì)能

會(huì)帶來(lái)非對(duì)角元。

定態(tài)薛定諤方程的數(shù)值解

這里我用 python 把一維定態(tài)薛定諤方程的數(shù)值解封裝成一個(gè)類，后面研究不同勢(shì)能下的薛定諤方程比較方便。

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):'''構(gòu)造二階微分算子：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):'''解哈密頓矩陣的本征值，本征向量；并對(duì)本征向量排序'''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$')

諧振子勢(shì)

# 定義諧振子勢(shì) def harmonic_potential(x, k=100):return 0.5 * k * x**2# 創(chuàng)建諧振子勢(shì)下的薛定諤方程 schro_harmonic = Schrodinger(harmonic_potential)

從上面例子可以看到，封裝的比較完整，對(duì)任意 1 維勢(shì)調(diào)用很簡(jiǎn)單。先來(lái)可視化諧振子勢(shì)能，

schro_harmonic.plot_potential()

schro_harmonic.check_eigen(n=1)

再用上面這條命令檢查一下薛定諤方程的解是否準(zhǔn)確，具體來(lái)說(shuō)就是本征方程

是否滿足。

還可以看看粒子在諧振子勢(shì)阱中的分布概率密度，

# 這里隨便選了一個(gè)能級(jí) n = 9 schro_harmonic.plot_density(n=9)

如果親自嘗試一下，你會(huì)發(fā)現(xiàn)在上面這個(gè)諧振子勢(shì)下，數(shù)值解的能級(jí)與解析解非常接近

對(duì)比解析解，

其中

, 。

解析解中，n=0 時(shí)，

。 n = 1, 2, 3... 時(shí)， 是 的 3 倍，5倍，7倍... 。

Woods Saxon 勢(shì)能

在核物理領(lǐng)域，原子核中一大團(tuán)核子所產(chǎn)生的勢(shì)能接近于 Woods Saxon 函數(shù)形式。這里看看Woods Saxon 勢(shì)阱中一個(gè)核子的能級(jí)分布。勢(shì)阱函數(shù)形式為，

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

用這個(gè)勢(shì)阱構(gòu)造薛定諤方程，

ws_schro = Schrodinger(woods_saxon_potential)

先畫(huà)一下勢(shì)阱的樣子，

ws_schro.plot_potential()

再畫(huà)一下 Woods Saxon 勢(shì)阱中核子的波函數(shù)和能級(jí)，

對(duì)于基態(tài) n=0

ws_schro.plot_density(n=0)

第 n=9 激發(fā)態(tài)

與諧振子勢(shì)阱的結(jié)果有相當(dāng)大差別。

雙勢(shì)阱

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)

雙勢(shì)阱中前幾個(gè)能級(jí)下粒子的概率密度分布，

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

雙勢(shì)阱中粒子概率密度分布的隨時(shí)間演化

下面這個(gè)問(wèn)題考慮一個(gè)粒子被捕獲在上例所示的有限深方勢(shì)阱中，初態(tài)為基態(tài)與第一激發(fā)態(tài)的疊加態(tài)，觀察粒子的概率密度分布隨時(shí)間的演化。初態(tài)為，

這里畫(huà)圖看一下基態(tài)

，第一激發(fā)態(tài) 和兩者疊加態(tài) 的波函數(shù)，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)

如下圖所示，初始時(shí)刻基態(tài)與第一激發(fā)態(tài)在右邊勢(shì)阱處相消，導(dǎo)致疊加態(tài)的波函數(shù)在左邊勢(shì)阱處有峰值結(jié)構(gòu)。后面演示此峰如何隨時(shí)間在兩個(gè)勢(shì)阱間振蕩。

疊加態(tài)波函數(shù)的時(shí)間演化直接用時(shí)間演化算符，

def psit(t, hbar=1):'''基態(tài)與第一激發(fā)態(tài)的疊加態(tài)波函數(shù)，隨時(shí)演化'''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

列向量表示所有離散空間點(diǎn)上的波函數(shù)值。用 表示給定 點(diǎn)上的波函數(shù)值。波函數(shù)的平方表示概率密度，對(duì)于給定的離散時(shí)空點(diǎn)

當(dāng)

函數(shù)項(xiàng)的值從 1 變到負(fù) 1 ， 點(diǎn)的概率密度從

變到

下面是概率密度在兩個(gè)勢(shì)阱間震蕩間振蕩的動(dòng)畫(huà)，

知乎視頻?www.zhihu.com# 動(dòng)畫(huà)代碼 %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,# 畫(huà)縮放了的雙勢(shì)阱 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()

總結(jié)：

使用微分的有限差分近似可以將波函數(shù)表示為向量，將微分算子化為矩陣，將定態(tài)薛定諤方程的求解化為哈密頓矩陣的本征值，本征向量求解問(wèn)題。實(shí)現(xiàn)了 python 版本的一維薛定諤方程數(shù)值解的封裝，方便對(duì)自定義的勢(shì)阱計(jì)算能級(jí)與波函數(shù)。

參考文獻(xiàn)：

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

【2】Teaching Quantum Mechanics with MATLAB

注：畫(huà)圖用到了 matplotlib 庫(kù)，要得到本文畫(huà)圖風(fēng)格，需要安裝 SciencePlots 庫(kù)。

pip install SciencePlots

如果出錯(cuò)，開(kāi)啟 no-latex 選項(xiàng)。

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

總結(jié)

以上是生活随笔為你收集整理的用diag直接使用错误_用python学量子力学（1）的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。