跳转至

谐振子数值解

我们知道,对于做谐振的粒子,它的薛定谔方程为:

\[ -\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2}+\frac{1}{2}kx^2=E \]

通过解这个薛定谔方程我们可以得到波函数的解析解,进而得到能量满足:

\[ E_n=(n+\frac{1}{2})\hbar\omega \]

现在,我们通过一个解析解的方式来尝试得出不同能级的能量值:

我们知道,波函数图像中的哈密顿算子为:

\[ \hat{H}=\hat{T}+\hat{V}=-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}+\frac{1}{2}kx^2 \]

换到矩阵图像中,我们分别对势能算子和动能算子展开讨论:

首先要选定一组基,这里我们采用坐标表象的基,也就是坐标算符的本征函数作为基:

注意到根据\(\delta\)函数的选择性有:

\[ f(x_0)=\int_{-\infty}^{+\infty} f(x)\delta(x-x_0)dx \]

现在让这个\(x_0\)跑动起来:

\[ f(t)=\int_{-\infty}^{+\infty} f(x)\delta(x-t)dx \]

这样就得到了一个积分变换,这是一个无穷维度的求和,可以写成黎曼和的形式(在一个有限的大区间\(\Delta L\)内):

\[ f(t)=\lim_{N \to \infty}\sum_{i=1}^N\frac{f(x_i)\Delta L}{N}\delta(t-x_i) \]

所以我们选定的基矢量就是:

\[ \left( |\delta(x-x_1)\rangle,|\delta(x-x_2)\rangle, \ldots ,|\delta(x-x_N)\rangle \right) \]

但是,这个\(\delta\)函数不是一个初等函数,我们不能很好的表示它,所以通常会用一些辅助函数来逼近,例如:

\[ \begin{aligned} f(x)=\frac{\sin x}{x} \\ y=\frac{1}{\sqrt{2\pi}}e^{\frac{x^2}{2}} \end{aligned} \]

势能算子:

\[ V_{ij}=\langle\psi_i|\frac{1}{2}kx^2|\psi_j\rangle=\int \delta(x-x_i) \frac{1}{2}kx^2 \delta(x-x_j)dx= \begin{cases} \frac{1}{2}kx_i^2, &i=j \\0 , & i\neq j \end{cases} \]

所以势能算子是对角阵:

\[ \hat{V}= \begin{pmatrix} &\frac{1}{2}kx_1^2,&0, &\ldots ,&0\\ &0,&\frac{1}{2}kx_2^{2}, &\ldots ,&0\\ &\vdots, &\vdots ,&\ddots,&\vdots\\ &0,&0,&\cdots,&\frac{1}{2}kx_N^2 \end{pmatrix} \]

对动能算子,我们考虑用有限差分法去处理二阶导数:

考虑用差分代替微分:

\[ \frac{df(x)}{dx}\approx \frac{f(x+\Delta x)-f(x)}{\Delta x} \]

对二阶导数有:

\[ \frac{d^2f(x)}{dx^2}\approx \frac{f^{'}(x+\frac{\Delta x}{2})-f^{'}(x-\frac{\Delta x}{2})}{\Delta x}=\frac{f(x+\Delta x)+f(x-\Delta x)-2f(x)}{\Delta x^2} \]

所以说:

\[ \frac{d^2}{dx^2}|\psi(x)\rangle=\frac{|\psi(x+\Delta x)\rangle+|\psi(x-\Delta x)\rangle-2|\psi(x)\rangle}{\Delta x^2} \]

如何定义我们的基矢量,考虑一个充分大的区间\([-L,L]\),这个充分大要到什么程度呢,即波函数基本衰减到0,这样我们能够用这组基尽可能准确的描述波函数的信息:

\[ \Psi(\pm L)\approx 0 \]

将区间划分成为\(N\)等份,取包括端点内的\(N+1\)个点:

\[ x_i=x_0+i \frac{2L}{N},\quad i=0,1,2, \ldots ,N \]

基矢量为:

\[ (|\delta(x-x_1)\rangle,|\delta(x-x_2)\rangle, \ldots ,|\delta(x-x_N)\rangle) \]

考虑两个相同的基矢量在算符作用下的内积:

\[ \langle\psi_i(x)|\frac{d^2}{dx^2}|\psi_i(x)\rangle=-\frac{2}{\Delta x^2}\langle \psi_i|\psi_i\rangle=-\frac{2}{\Delta x^2} \]

对于两个不同的基矢量,要进行分类讨论:

如果说:\(x_j=x_i\pm \Delta x\),即\(j=i\pm 1\):

\[ \langle\psi_j(x)|\frac{d^2}{dx^2}|\psi_i(x)\rangle=\langle\psi_j(x)|\frac{|\psi(x+\Delta x)\rangle+|\psi(x-\Delta x)\rangle-2|\psi(x)\rangle}{\Delta x^2}= \frac{1}{\Delta x^2} \]

其余的矩阵元由于不相邻,内积都为0.

所以动量算子的矩阵为:

\[ \hat{T}=-\frac{\hbar^2}{2m} \begin{pmatrix} &-\frac{2}{\Delta x^2},&\frac{1}{\Delta x^2}, &\ldots ,&0,&0\\ &\frac{1}{\Delta x^2},&-\frac{2}{\Delta x^2},&\frac{1}{\Delta x^2}, &\ldots ,&0\\ &\vdots, &\vdots ,&\ddots,&\ddots,&\vdots\\ &0,&0,&\cdots,&\frac{1}{\Delta x^2},&-\frac{2}{\Delta x^2} \end{pmatrix} \]

总的哈密顿算符就由这两个算符相加得到一个新矩阵,下一步就是把新矩阵对角化,对角线上的值就是矩阵的特征值,这个特征值就是我们要求的谐振子能量的数值解,具体可以通过Python的numpy库来实现:

代码如下:

import numpy as np

#定义维度数

N=1000

#把常数k,m,hbar,w全部定义成1

#定义区间和划分区间

a=-10
b=10
x=[]
for i in range(N):
    x.append(a+(b-a)/N*(i+1))

# 定义势能矩阵

V=np.zeros((N,N))

for i in range(N):
    for j in range(N):
        if i==j:
            V[i][j]=0.5*x[i]**2

#定义动能矩阵

T=np.zeros((N,N))

dx=(b-a)/N

for i in range(N):
    for j in range(N):
        if i==j:
            T[i][j]=-2/dx**2
        elif j==i+1 or j==i-1:
            T[i][j]=1/dx**2

#哈密顿算子

H=V+(-1/2)*T

# 求解特征值

eigenvalues,eigenvectors= np.linalg.eig(H)

eigenvalues.sort()

print('前几个特征值:')

for i in range(5):
    print(f'{eigenvalues[i]:}')

运行结果如下:

前几个特征值:
0.4999874996876701
1.4999374971917387
2.499837489063911
3.499687471559114
4.499487440926194

与解析解符合相当完美,精度足够.