Finite Difference Approximation
The second derivative is approximated using central differences:
\[
\frac{d^2\psi}{dx^2}\Big|_{x_i} \approx
\frac{\psi_{i+1} - 2\psi_i + \psi_{i-1}}{(\Delta x)^2}
\]
where \( \Delta x = L/(N - 1) \) is the grid spacing.
Hamiltonian in Matrix Form
Discretizing the spatial domain into \( N \) points, the wavefunction becomes a column vector:
\[
\vec{\psi} =
\begin{bmatrix}
\psi_1 \\[4pt]
\psi_2 \\[2pt]
\vdots \\[2pt]
\psi_N
\end{bmatrix}
\]
The kinetic energy operator using finite differences is represented by the tridiagonal matrix \( T \):
\[\small
T = \frac{1}{(\Delta x)^2}
\begin{bmatrix}
-2 & 1 & 0 & \cdots & 0 \\
1 & -2 & 1 & \cdots & 0 \\
0 & 1 & -2 & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & 1 & -2
\end{bmatrix}
\]
The potential energy operator becomes a diagonal matrix:
\[\small
V =
\begin{bmatrix}
V_1 & 0 & \cdots & 0 \\
0 & V_2 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & V_N
\end{bmatrix}
\]
The Hamiltonian matrix is then:
\[
H = -\frac{\hbar^2}{2m} T + V
\]
The time-independent Schrödinger equation in matrix form is:
\[
H \vec{\psi} = E \vec{\psi}
\]
or explicitly, component-wise:
\[\tiny
\begin{bmatrix}
H_{11} & H_{12} & 0 & \cdots & 0 \\
H_{21} & H_{22} & H_{23} & \cdots & 0 \\
0 & H_{32} & H_{33} & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & H_{N-1,N} \\
0 & 0 & 0 & H_{N,N-1} & H_{NN}
\end{bmatrix}
\begin{bmatrix}
\psi_1 \\ \psi_2 \\ \vdots \\ \psi_N
\end{bmatrix}
=
E
\begin{bmatrix}
\psi_1 \\ \psi_2 \\ \vdots \\ \psi_N
\end{bmatrix}
\]
where
\[\small
H_{ii} = \frac{\hbar^2}{2m(\Delta x)^2} + V_i, \quad
H_{i,i\pm1} = -\frac{\hbar^2}{2m(\Delta x)^2},
\]
and all other elements are zero.
Solving this eigenvalue problem yields the allowed energies \( E \) and the corresponding normalized wavefunctions \( \vec{\psi} \).
Expectation Values
Expectation values are computed numerically as:
\[
\langle x \rangle = \int x |\psi(x)|^2 \, dx
\]
\[
\langle x^2 \rangle = \int x^2 |\psi(x)|^2 \, dx
\]
Algorithm
- Define the potential \( V(x) \) on a discrete grid of \( N \) points over domain \( L \).
- Construct the kinetic energy matrix \( T \) using finite differences (tridiagonal matrix).
- Construct the Hamiltonian matrix: \( H = -\frac{\hbar^2}{2m}T + \text{diag}(V) \).
- Solve the eigenvalue problem \( H\vec{\psi} = E\vec{\psi} \) using
numpy.linalg.eigh.
- Normalize the wavefunctions.
- Plot probability densities \( |\psi_n(x)|^2 \).
- Compute expectation values \( \langle x \rangle \) and \( \langle x^2 \rangle \) using the trapezoidal rule.