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

## 一维谐振子

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 ^ ∣ ψ > = E ∣ ψ > \hat H\left|\psi\right> = E\left|\psi\right>

## 有限差分法

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)
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 ′ ′ ( 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}

ψ ( x ) \psi(x) x ∈ [ − r , r ] x\in[-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

− ℏ 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

[ 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]

## QR 算法

QR 算法是一种常见的特征值算法。它利用了矩阵的 QR 分解，即将矩阵 A A 分解为一个正交矩阵 Q Q 和一个上三角矩阵 R R 的乘积：
A = Q R A=QR

A 0 : = A A_0:=A ，对 k = 0 , 1 , 2 , ⋯ k=0, 1, 2, \cdots
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}

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}

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 k A_k 收敛时， S k S_k 的列向量即为属于相应特征值的特征向量。

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}
A k A_k 收敛时，就同时求得 A A 的特征值和特征向量。

## Code

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


# 取 [-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]


## 与解析解的对比

E n = ( n + 1 2 ) ℏ ω , n = 0 , 1 , 2 , ⋯ E_n = \left(n+\frac12\right)\hbar\omega,\qquad n=0, 1, 2,\cdots

ψ 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)

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)

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}

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}

## 参考文献

UP更新不错过~
• 3
点赞
• 16
收藏
• 打赏
• 8
评论
12-05
01-11 247
03-10 1382
10-27 2945
03-09 632
07-22 680
12-18 5492
11-10 3630
01-10 6459
01-02

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

• 非常没帮助
• 没帮助
• 一般
• 有帮助
• 非常有帮助

¥2 ¥4 ¥6 ¥10 ¥20

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