Solving the Schrödinger Equation by exact diagonalization

I am solving the Schrödinger equation via finite difference, via the substitution

where we are assuming $$V_1 = V_N = \infty$$. I solved this using Mathematica for the case that $$V(x) = 0$$ and get the correct eigenvalues (in units with $$\hbar = m = 1$$, measuring length in units of the width of the well $$a$$, the eigenvalues are simply $$n^2 \pi^2$$). I find the eigenvectors and plot them via

d2dx[n_] :=
SparseArray[{Band[{1, 1}] -> -2, Band[{1, 2}] -> 1,
Band[{2, 1}] -> 1}, {n, n}]/(1/n^2);
evecs[n_, m_] := Last@Eigenvectors[N[-d2dx[n]], m]
ListPlot[evecs[100, 1], Joined -> True]
ListPlot[evecs[100, 2], Joined -> True]
ListPlot[evecs[100, 3], Joined -> True]


But the results look like this:

The profiles are correct but of course it shouldn't be oscillating in this way. Every other point is reflected. I am not sure if there is an issue in how I am solving for the eigenvectors or if this is a mathematical issue. Does anyone know why am I seeing these oscillations and how do I get the correct eigenvectors?

• How do you know those aren't correct eigenvectors? Did you try $\left\langle \psi \middle| H \middle| \psi \right\rangle$ to see if it returns the associated eigenvalue? Generally, eigenvectors are only defined up to a constant phase, so I don't see that these are necessarily wrong. May 21 '18 at 18:03
• As @JasonB. notes, commonly for things like this you need to correct with a known phase. I.e. you need to make sure the ground-state is node-less and then apply the same phase correction to everything else. The process you're using is akin to the idea behind quadrature DVRs, where you most certainly need to do that. Oct 22 '18 at 8:03

Eigenvalues and Eigenvectors return the eigenvalues sorted from highest to lowest. From the documentation:

Eigenvectors with numeric eigenvalues are sorted in order of decreasing absolute value of their eigenvalues.

This means that the last three eigenvectors returned by Eigenvectors will correspond to the lowest-energy states:

ListPlot[{evecs[100, 98], evecs[100, 99], evecs[100, 100]}, Joined -> True]


The eigenvectors you found were valid solutions to the discretized problem. However, given that your "grid spacing" was comparable to the "wavelength" of the particle for these functions, it is doubtful that they were physically meaningful.

Addendum: Your original problem is equivalent to the normal modes of a beaded string, for which the eigenvalues are known to be $$\omega^2 \propto 4 \sin^2 \left( \frac{r \pi}{2 (n+1)} \right), \qquad r = 1, 2, \dots, n$$ where $n$ is the number of "beads" on the string. If you're familiar with this system, the patterns you found corresponding to the highest eigenvalues (i.e. highest normal frequencies) might make a bit more sense. For more information, see the link above; Thornton & Marion's Classical Dynamics of Particles & Systems has a more in-depth discussion of the system.

• Are these wavefunctions normalized? May 26 '18 at 9:44
• @AnkurMandal: From the documentation: "For approximate numerical matrices m, the eigenvectors are normalized." Since the code invokes N before finding the eigenvectors, they're therefore normalized. That said, "normalized" in this context means that the sum of the squares of the discretized entries is 1; I think this is equivalent to saying that $\int |\psi^2| dx = 1$, but it depends on the normalization factors and on the way the code is written. To be sure you'd have to approximate the integral by a Riemann sum and compare to the normalization condition above. May 26 '18 at 13:33