一维谐振子定态 Schrödinger 方程的数值解法

本文链接个人站 | 简书 | CSDN
版权声明:除特别声明外,本博客文章均采用 BY-NC-SA 许可协议。转载请注明出处。

前几天整理电脑的时候发现了本科上量子力学讨论班时做的一个 Slide,觉得挺有意思的。花了点时间整理成这篇博客。

一维谐振子

一个质量为 m m m 的粒子,在一维势场 V ( x ) = 1 2 m ω 2 x 2 V(x) = \dfrac12m\omega^2x^2 V(x)=21mω2x2 中运动。其哈密顿算符为
H ^ = p ^ 2 2 m + 1 2 m ω 2 x ^ 2 \hat H = \frac{\hat p^2}{2m} + \frac12m\omega^2\hat x^2 H^=2mp^2+21mω2x^2
其中 x ^ \hat x x^ 为位置算符, p ^ = − i ℏ d d x \hat p = -i\hbar\dfrac{\mathrm d}{\mathrm dx} p^=idxd 为动量算符。我们需要求解该体系的定态 Schrödinger 方程:
H ^ ∣ ψ > = E ∣ ψ > \hat H\left|\psi\right> = E\left|\psi\right> H^ψ=Eψ

一维谐振子是除了氢原子之外,为数不多的可以解析求解的体系。那么我们为什么要费劲求它的数值解呢?正因为绝大多数的量子体系都无法解析求解,数值方法才显得尤为重要。

有限差分法

回忆一下泰勒公式
f ( a + h ) = f ( a ) + f ′ ( a ) 1 ! h + f ′ ′ ( a ) 2 ! h 2 + o ( h 3 ) f(a+h) = f(a) + \frac{f'(a)}{1!}h+\frac{f''(a)}{2!}h^2+o(h^3) f(a+h)=f(a)+1!f(a)h+2!f(a)h2+o(h3)
h = − h h=-h h=h,有
f ( a − h ) = f ( a ) − f ′ ( a ) 1 ! h + f ′ ′ ( a ) 2 ! h 2 + o ( h 3 ) f(a-h) = f(a) - \frac{f'(a)}{1!}h+\frac{f''(a)}{2!}h^2+o(h^3) f(ah)=f(a)1!f(a)h+2!f(a)h2+o(h3)
两式相加,可得
f ′ ′ ( a ) = f ( a − h ) + f ( a + h ) − 2 f ( a ) h 2 + o ( h 3 ) ≈ f ( a − h ) + f ( a + h ) − 2 f ( a ) h 2 \begin{aligned} f''(a) &= \frac{f(a-h)+f(a+h)-2f(a)}{h^2} + o(h^3)\\ &\approx \frac{f(a-h)+f(a+h)-2f(a)}{h^2} \end{aligned} f(a)=h2f(ah)+f(a+h)2f(a)+o(h3)h2f(ah)+f(a+h)2f(a)

ψ ( x ) \psi(x) ψ(x) x ∈ [ − r , r ] x\in[-r, r] x[r,r] 区间离散化为
ϕ i ≡ ψ ( x i ) = ψ ( i Δ x − r ) , i = 0 , 1 , 2 , ⋯   , N \phi_i \equiv \psi(x_i) = \psi(i\Delta x - r),\quad i=0, 1, 2, \cdots, N ϕiψ(xi)=ψ(iΔxr),i=0,1,2,,N
其中 N = 2 r / Δ x N=2r/\Delta x N=2r/Δx,则 Schrödinger 方程差分化为
− ℏ 2 2 m ϕ i − 1 + ϕ i + 1 − 2 ϕ i Δ x 2 + 1 2 m ω 2 x i 2 ϕ i = E ϕ i -\frac{\hbar^2}{2m}\frac{\phi_{i-1}+\phi_{i+1}-2\phi_i}{\Delta x^2}+\frac12m\omega^2x_i^2\phi_i = E\phi_i 2m2Δx2ϕi1+ϕi+12ϕi+21mω2xi2ϕi=Eϕi
在这里,我们假设 x < − r x<-r x<r x > r x>r x>r ψ ( x ) → 0 \psi(x)\to 0 ψ(x)0。这对能量较低的态是成立的。

将差分方程写成矩阵的形式为
[ m ω 2 x 0 2 2 + ℏ 2 m Δ x 2 − ℏ 2 2 m Δ x 2 0 ⋯ 0 − ℏ 2 2 m Δ x 2 m ω 2 x 1 2 2 + ℏ 2 m Δ x 2 − ℏ 2 2 m Δ x 2 ⋯ 0 0 − ℏ 2 2 m Δ x 2 m ω 2 x 2 2 2 + ℏ 2 m Δ x 2 ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ m ω 2 x N 2 2 + ℏ 2 m Δ x 2 ] [ ϕ 0 ϕ 1 ϕ 2 ⋮ ϕ N ] = E [ ϕ 0 ϕ 1 ϕ 2 ⋮ ϕ N ] \left[ \begin{matrix} \frac{m\omega^2x_0^2}2+\frac{\hbar^2}{m\Delta x^2} & -\frac{\hbar^2}{2m\Delta x^2} & 0 & \cdots & 0\\ -\frac{\hbar^2}{2m\Delta x^2} & \frac{m\omega^2x_1^2}2+\frac{\hbar^2}{m\Delta x^2} & -\frac{\hbar^2}{2m\Delta x^2} & \cdots & 0\\ 0 & -\frac{\hbar^2}{2m\Delta x^2} & \frac{m\omega^2x_2^2}2+\frac{\hbar^2}{m\Delta x^2} & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \frac{m\omega^2x_N^2}2+\frac{\hbar^2}{m\Delta x^2} \end{matrix} \right] \left[ \begin{matrix} \phi_0\\ \phi_1\\ \phi_2\\ \vdots\\ \phi_N \end{matrix} \right]= E \left[ \begin{matrix} \phi_0\\ \phi_1\\ \phi_2\\ \vdots\\ \phi_N \end{matrix} \right] 2mω2x02+mΔx222mΔx22002mΔx222mω2x12+mΔx222mΔx22002mΔx222mω2x22+mΔx2200002mω2xN2+mΔx22ϕ0ϕ1ϕ2ϕN=Eϕ0ϕ1ϕ2ϕN

这样一来,问题就转化为求差分矩阵的特征值和特征向量。

QR 算法

QR 算法是一种常见的特征值算法。它利用了矩阵的 QR 分解,即将矩阵 A A A 分解为一个正交矩阵 Q Q Q 和一个上三角矩阵 R R R 的乘积:
A = Q R A=QR A=QR
为什么可以这么分解呢?我们回忆一下 Gram-Schmidt 正交化,将矩阵 A A A 的列向量看作一组基,则可以通过一系列初等列变换获得一组标准正交基。反过来看,对这组标准正交基所组成的矩阵 Q Q Q 作初等列变换也可以得到矩阵 A A A。我们知道,对矩阵作初等列变换相当于右乘初等矩阵,且由于 Gram-Schmidt 正交化不涉及列交换,这里用到的初等矩阵均为上三角矩阵。因此,矩阵 A A A 可以由正交矩阵 Q Q Q 右乘一个上三角矩阵 R R R 得到。在实际应用中,除了 Gram-Schmidt 正交化,还可以用 Householder 变换Givens 旋转 等方法实现 QR 分解。

那么如何利用 QR 分解求解矩阵的特征值呢?

A 0 : = A A_0:=A A0:=A,对 k = 0 , 1 , 2 , ⋯ k=0, 1, 2, \cdots k=0,1,2,
A k = Q k R k A k + 1 : = R k Q k \begin{aligned} A_k &= Q_kR_k\\ A_{k+1} &:= R_kQ_k \end{aligned} AkAk+1=QkRk:=RkQk
即在每一步,对 A k A_k Ak 进行 QR 分解,再由分解后得到对 Q k Q_k Qk R k R_k Rk 计算 A k + 1 A_{k+1} Ak+1,如此迭代。

注意到 Q k Q_k Qk 是正交矩阵,有 Q k − 1 = Q k ⊤ Q_k^{-1}=Q_k^\top Qk1=Qk,故
A k + 1 = R k Q k = Q k − 1 Q k R k Q k = Q k − 1 A k Q k = Q k ⊤ A k Q k \begin{aligned} A_{k+1} &= R_kQ_k\\ &= Q_k^{-1} Q_kR_kQ_k\\ &= Q_k^{-1}A_kQ_k\\ &= Q_k^\top A_kQ_k \end{aligned} Ak+1=RkQk=Qk1QkRkQk=Qk1AkQk=QkAkQk
也就是说, A k + 1 A_{k+1} Ak+1 相似于 A k A_k Ak。根据递推关系, A 0 , A 1 , ⋯   , A k , ⋯ A_0, A_1, \cdots, A_k, \cdots A0,A1,,Ak, 全都是相似的,这意味着所有的 A k A_k Ak 都有相同的特征值。在一定条件下, A k A_k Ak 会收敛为一个三角矩阵,特征值为其主对角元。

特别地,如果 A A A 是一个实对称正定矩阵(我们所要求解的差分矩阵刚好是这种情况), A k A_k Ak 将会收敛为一个对角矩阵 Λ = d i a g { λ 0 , λ 1 , ⋯   , λ N } \Lambda=diag\{\lambda_0,\lambda_1,\cdots,\lambda_N\} Λ=diag{λ0,λ1,,λN},且 [ λ 0 , λ 1 , ⋯   , λ N ] [\lambda_0,\lambda_1,\cdots,\lambda_N] [λ0,λ1,,λN] 依次递减。考虑到
A = A 0 = Q 0 A 1 Q 0 ⊤ = ( Q 0 Q 1 ⋯ Q k − 1 ) A k ( Q k − 1 ⊤ ⋯ Q 1 ⊤ Q 0 ⊤ ) = ( Q 0 Q 1 ⋯ Q k − 1 ) A k ( Q 0 Q 1 ⋯ Q k − 1 ) ⊤ = S k A k S k ⊤ \begin{aligned} A=A_0 &= Q_0A_1Q_0^\top\\ &=(Q_0 Q_1\cdots Q_{k-1})A_k(Q_{k-1}^\top\cdots Q_1^\top Q_0^\top)\\ &=(Q_0 Q_1\cdots Q_{k-1})A_k(Q_0Q_1\cdots Q_{k-1})^\top\\ &=S_kA_kS_k^\top \end{aligned} A=A0=Q0A1Q0=(Q0Q1Qk1)Ak(Qk1Q1Q0)=(Q0Q1Qk1)Ak(Q0Q1Qk1)=SkAkSk
A k A_k Ak 收敛时, S k S_k Sk 的列向量即为属于相应特征值的特征向量。

综上,对于实对称正定矩阵 A A A,我们令 A 0 = A , S 0 = I A_0=A, S_0=I A0=A,S0=I,对 k = 0 , 1 , 2 , ⋯ k=0, 1, 2, \cdots k=0,1,2,
A k = Q k R k A k + 1 : = R k Q k S k + 1 : = S k Q k \begin{aligned} A_k &= Q_kR_k\\ A_{k+1} &:= R_kQ_k\\ S_{k+1} &:= S_kQ_k \end{aligned} AkAk+1Sk+1=QkRk:=RkQk:=SkQk
A k A_k Ak 收敛时,就同时求得 A A A 的特征值和特征向量。

Code

在 Python 下可以直接调用 numpy.linalg.qr 作 QR 分解:

import numpy as np

def qr_eig(A, iters=200, tol=1e-6):
    """
    使用 QR 算法求解实对称矩阵的特征值和特征向量
    
    Parameters
    ----------
    A : Array-like
        二维数组,表示待求解的实对称矩阵
    iters : int
        最大迭代次数
    tol : float
        提前退出循环的判断标准
    
    Returns
    -------
    (numpy.ndarray, numpy.ndarray)
        一维数组表示的特征值,二维数组表示的特征向量
    """
    A = np.asarray(A)
    S = np.eye(A.shape[0])
    for _ in range(iters):
        Q, R = np.linalg.qr(A)
        newA = np.dot(R, Q)
        newS = np.dot(S, Q)
        if np.abs(max(np.diag(newA)) - max(np.diag(A))) < tol:
            break
        A = newA
        S = newS
    return np.diag(A), S

接下来就可以求解一维谐振子了。为了方便起见,我们令 ℏ = m = ω = 1 \hbar=m=\omega=1 =m=ω=1

# 取 [-5, 5] 的区间,离散化为 N 个点
N = 101
r = 5
dx = 2 * r / (N - 1)
x = np.linspace(-r, r, N, endpoint=True)

# 初始化差分矩阵
A = np.diag(0.5*x**2 + 1/dx**2)
for i in range(N -1 ):
    A[i][i+1] = -0.5/dx**2
    A[i+1][i] = -0.5/dx**2

# 求解差分矩阵的特征值和特征向量
Lambda, S = qr_eig(A)

# 结果展示
# 打印最低的 5 个能级
print(Lambda[-1:-8:-1])
# 画出能量最低的三个态的波函数
plt.plot(x, S[:,-1], label=f'n=0, E={Lambda[-1]:.4f}')
plt.plot(x, S[:,-2], label=f'n=1, E={Lambda[-2]:.4f}')
plt.plot(x, S[:,-3], label=f'n=2, E={Lambda[-3]:.4f}')
plt.xlabel('x')
plt.ylabel('psi(x)')
plt.legend()
>>> [0.4996873  1.49843574 2.49593067 3.49217016 4.48715598]

numerical.png

与解析解的对比

一维谐振子的本征能量为:
E n = ( n + 1 2 ) ℏ ω , n = 0 , 1 , 2 , ⋯ E_n = \left(n+\frac12\right)\hbar\omega,\qquad n=0, 1, 2,\cdots En=(n+21)ω,n=0,1,2,
对应的本征态为:
ψ n ( x ) = 1 2 n n ! ⋅ ( m ω π ℏ ) 1 / 4 ⋅ e − m ω x 2 / 2 ℏ ⋅ H n ( m ω ℏ x ) \psi_n(x) = \frac{1}{\sqrt{2^n n!}}\cdot\left(\frac{m\omega}{\pi\hbar}\right)^{1/4}\cdot\mathrm{e}^{-m\omega x^2/2\hbar}\cdot H_n\left(\sqrt{\frac{m\omega}{\hbar}}x\right) ψn(x)=2nn! 1(πmω)1/4emωx2/2Hn(mω x)
其中
H n ( z ) = ( − ) n e z 2 d n d z n ( e − z 2 ) H_n(z) = (-)^n\mathrm{e}^{z^2}\frac{\mathrm d^n}{\mathrm dz^n}\left(\mathrm{e}^{-z^2}\right) Hn(z)=()nez2dzndn(ez2)
为 Hermite 多项式。前三个 Hermite 多项式为:
H 0 ( z ) = 1 H 1 ( z ) = 2 z H 2 ( z ) = 4 z 2 − 2 \begin{aligned} H_0(z) &= 1\\ H_1(z) &= 2z\\ H_2(z) &= 4z^2 - 2 \end{aligned} H0(z)H1(z)H2(z)=1=2z=4z22
我们同样令 ℏ = m = ω = 1 \hbar=m=\omega=1 =m=ω=1,则
E 0 = 1 2 , ψ 0 ( x ) = 1 π 1 / 4 ⋅ e − x 2 / 2 E 1 = 3 2 , ψ 1 ( x ) = 1 2 ⋅ 1 π 1 / 4 ⋅ e − x 2 / 2 ⋅ 2 x E 2 = 5 2 , ψ 2 ( x ) = 1 2 ⋅ 1 π 1 / 4 ⋅ e − x 2 / 2 ⋅ ( 2 x 2 − 1 ) \begin{aligned} E_0&=\frac12, \qquad \psi_0(x) = \frac{1}{\pi^{1/4}}\cdot\mathrm{e}^{-x^2/2}\\ E_1&=\frac32, \qquad \psi_1(x) = \frac{1}{\sqrt 2}\cdot \frac{1}{\pi^{1/4}}\cdot\mathrm{e}^{-x^2/2}\cdot 2x\\ E_2&=\frac52, \qquad \psi_2(x) = \frac{1}{\sqrt 2}\cdot \frac{1}{\pi^{1/4}}\cdot\mathrm{e}^{-x^2/2}\cdot (2x^2-1) \end{aligned} E0E1E2=21,ψ0(x)=π1/41ex2/2=23,ψ1(x)=2 1π1/41ex2/22x=25,ψ2(x)=2 1π1/41ex2/2(2x21)
画出来看看

analytical.png

八九不离十吧。低能级的误差主要来自截断误差和舍入误差。此外,高能级需要有更大的 r r r 来保证 x < − r x<-r x<r x > r x>r x>r ψ ( x ) → 0 \psi(x)\to 0 ψ(x)0 的假设,因此能级越高,误差越大。

参考文献

  1. Quantum harmonic oscillator - Wikipedia
  2. Finite difference method - Wikipedia
  3. QR algorithm - Wikipedia
  4. Notes on orthogonal bases and the workings of the QR algorithm
  5. QR decomposition - Wikipedia
  6. Gram–Schmidt process - Wikipedia
微信扫码订阅
UP更新不错过~
关注
  • 3
    点赞
  • 16
    收藏
  • 打赏
    打赏
  • 8
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:数字20 设计师:CSDN官方博客 返回首页
评论 8

打赏作者

虚胖一场

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值