Finite element method for the 'Particle-In-a-Box' problem in quantum mechanics

Question

(Apologies in advance for the lengthy question, but it really is needed for a precise description of what I've done!)

In suitable units, the 'Particle-in-a-box' problem is described by the following time-independant Schrodinger equation and accompanying boundary conditions:

$-\frac{1}{2}\frac{d^2\psi}{dx^2} = E\psi; x \in [0, 1]$

$\psi(0) = \psi(1) = 0$.

I have attempted to solve the equation numerically (as a toy problem) by replacing $\frac{d^x}{dx^2}$ with a finite difference, resulting in the following recurrence relation:

$-a\psi_{n - 1} + 2a\psi_{n} -a\psi_{n + 1} = E \psi_{n};a = \frac{1}{2\Delta{x}^2}$

which is equivalent to the following system of N linear equations in $\psi_{n}$ (when the domain of x is divided into N points $x_n$):

$-a\psi_{0} + 2a\psi_{1} - a\psi_{2} = E \psi_{1}$

$-a\psi_{1} + 2a\psi_{2} - a\psi_{2} = E \psi_{2}$

$-a\psi_{2} + 2a\psi_{3} - a\psi_{4} = E \psi_{2}$

...

$-a\psi_{N - 2} + 2a\psi_{N - 1} - a\psi_{N} = E \psi_{N - 1}$

$-a\psi_{N - 1} + 2a\psi_{N} - a\psi_{N + 1} = E \psi_{N}$.

The system of equations is equivalent to the following matrix eigenvalue problem (e.g. when the section domain between the boundaries is divided into five points):

$$H|\psi> = E|\psi>$$

$$H = \begin{bmatrix} 2a & -a & 0 & 0 & 0 \\ -a & 2a & -a & 0 & 0 \\ 0 & -a & 2a & -a & 0 \\ 0 & 0 & -a & 2a & -a \\ 0 & 0 & 0 & -a & 2a \\ \end{bmatrix}$$

$$|\psi> = \begin{bmatrix} \psi_{1} \\ \psi_{2} \\ \psi_{3} \\ \psi_{4} \\ \psi_{5} \\ \end{bmatrix}$$

(Where the choice $\psi_{0} = \psi_{N + 1} = 0$ encodes the boundary conditions)

I have written the following python code to solve the problem in the case where N = 100:

import numpy as np

# input parameters
steps = 100

# internal parameters
dx = 1.0 / float(steps - 1)
a = 1 / float(2 * dx * dx)

# prepare hamiltonian matrix
hamiltonian = np.zeros((steps, steps))

for i in range(steps):
  for j in range(steps):
    if i == j:
      hamiltonian[i, j] = 2 * a
    elif i == j + 1 or j == i + 1:
      hamiltonian[i, j] = -a
# compute and unpack results
values, vectors = np.linalg.eig(hamiltonian)

The computed eigenvalues are out by many orders of magnitude (e.g. expected $E_{1} \approx E_{1}^{analytical} = \frac{\pi^2}{2} = 4.9348$ (to 5sf), but got $E_{1} = 10943$ (to 5sf) and the eigenvectors, when plotted against x appeared smoothly modulated and rapidly oscillating (nothing like the expected $\sqrt(2) * sin(\pi x)$). I tried again from scratch again and again but I cannot for the life of me see what I've been doing wrong. Can anyone help me? Thank-you in advance!

Answer 1

I ran the following Mathematica code:
n = 30;
H = IdentityMatrix[n]*2 - DiagonalMatrix[ConstantArray[1, n - 1], -1] - DiagonalMatrix[ConstantArray[1, n - 1], 1];
ListPlot[N[Eigenvectors[H]][[-5 ;; -1]], Joined -> True]
This creates the Hamiltonian (without the constant factor $a$) and plots the first eigenvectors
As you can see, is works as expected, so the error should be somewhere in you code. (The creation of $H$ seems fine, so I expect it has to do something with the linear algebra part in the last line.) Or, as mentioned as a comment, it has to do something with the expected energies of the system in the theoretical result.