華中師范大學 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"]):

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。