python矩阵方法解三维谐振子的薛定谔方程
时间: 2023-10-24 16:39:42 浏览: 75
三维谐振子的薛定谔方程为:
$$
-\frac{\hbar^2}{2m}\nabla^2\psi+\frac{1}{2}m\omega^2(x^2+y^2+z^2)\psi=E\psi
$$
其中 $m$ 为质量,$\omega$ 为振动频率,$E$ 为能量,$\psi$ 为波函数。
为了求解该方程,可以将其转化为矩阵形式,然后利用矩阵方法进行求解。
首先,将 $x,y,z$ 分别离散化为 $N$ 个点,令 $\Delta x = \Delta y = \Delta z = a$,则有:
$$
x_i = -\frac{aN}{2}+ia,\quad y_j = -\frac{aN}{2}+ja,\quad z_k = -\frac{aN}{2}+ka
$$
其中 $i,j,k$ 均为整数,$0\leq i,j,k<N$。
然后,将波函数离散化为 $N^3$ 个点,令 $\psi_{i,j,k}=\psi(x_i,y_j,z_k)$,则有:
$$
\nabla^2\psi_{i,j,k}=\frac{\psi_{i-1,j,k}+\psi_{i+1,j,k}+\psi_{i,j-1,k}+\psi_{i,j+1,k}+\psi_{i,j,k-1}+\psi_{i,j,k+1}-6\psi_{i,j,k}}{a^2}
$$
将上式代入原方程,得到:
$$
-\frac{\hbar^2}{2ma^2}\left(\psi_{i-1,j,k}+\psi_{i+1,j,k}+\psi_{i,j-1,k}+\psi_{i,j+1,k}+\psi_{i,j,k-1}+\psi_{i,j,k+1}-6\psi_{i,j,k}\right)+\frac{1}{2}m\omega^2(x_i^2+y_j^2+z_k^2)\psi_{i,j,k}=E\psi_{i,j,k}
$$
将上式化简,得到:
$$
-\frac{\hbar^2}{2ma^2}\left(\psi_{i-1,j,k}+\psi_{i+1,j,k}+\psi_{i,j-1,k}+\psi_{i,j+1,k}+\psi_{i,j,k-1}+\psi_{i,j,k+1}-6\psi_{i,j,k}\right)+\frac{1}{2}m\omega^2a^2\left(\frac{(i-N/2)^2+(j-N/2)^2+(k-N/2)^2}{4}\right)\psi_{i,j,k}=E\psi_{i,j,k}
$$
将 $\psi_{i,j,k}$ 放到一个 $N^3$ 维列向量 $\mathbf{\psi}$ 中,得到:
$$
H\mathbf{\psi}=E\mathbf{\psi}
$$
其中 $H$ 为 $N^3\times N^3$ 的矩阵,$\mathbf{\psi}$ 为 $N^3$ 维列向量,具体形式为:
$$
H = -\frac{\hbar^2}{2ma^2}\left(\mathbf{T}_x+\mathbf{T}_y+\mathbf{T}_z\right)+\frac{1}{2}m\omega^2a^2\mathbf{V}
$$
$$
\mathbf{T}_x=\begin{pmatrix}
-2 & 1 & 0 & \cdots & 0 \\
1 & -2 & 1 & \cdots & 0 \\
0 & 1 & -2 & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & -2 \\
\end{pmatrix}_{N\times N}
$$
$$
\mathbf{T}_y=\begin{pmatrix}
\mathbf{T}_x & \mathbf{0} & \cdots & \mathbf{0} \\
\mathbf{0} & \mathbf{T}_x & \cdots & \mathbf{0} \\
\vdots & \vdots & \ddots & \vdots \\
\mathbf{0} & \mathbf{0} & \cdots & \mathbf{T}_x \\
\end{pmatrix}_{N^2\times N^2}
$$
$$
\mathbf{T}_z=\begin{pmatrix}
\mathbf{T}_y & \mathbf{0} & \cdots & \mathbf{0} \\
\mathbf{0} & \mathbf{T}_y & \cdots & \mathbf{0} \\
\vdots & \vdots & \ddots & \vdots \\
\mathbf{0} & \mathbf{0} & \cdots & \mathbf{T}_y \\
\end{pmatrix}_{N^3\times N^3}
$$
$$
\mathbf{V}=\begin{pmatrix}
d_{1,1} & d_{1,2} & \cdots & d_{1,N} \\
d_{2,1} & d_{2,2} & \cdots & d_{2,N} \\
\vdots & \vdots & \ddots & \vdots \\
d_{N,1} & d_{N,2} & \cdots & d_{N,N} \\
\end{pmatrix}_{N\times N}
$$
其中 $\mathbf{0}$ 为全零矩阵,$d_{i,j}=\frac{(i-N/2)^2+(j-N/2)^2+(k-N/2)^2}{4}$。
最后,可以使用 numpy 中的 linalg.eigh 函数求解 $H\mathbf{\psi}=E\mathbf{\psi}$ 的本征值和本征向量,即可得到三维谐振子的能级和波函数。