How do I set up a numerical quantum particle in a box problem with an arbitrary potential?

Say I have the usual particle in a box problem, with a Hamiltonian $H=-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}+V(x)$, where $V$ is an arbitrary potential and the boundary conditions are such that $\psi(0)=0$ and $\psi(L)=0$. I know that for $V=0$ the eigensystem is described by eigenvalues

$$E_n=\frac{n^2 \pi^2 \hbar^2}{2m L^2}$$ and eigenvectors $$\langle x|\psi_n\rangle=\sqrt{\frac{2}{L}}\sin\left(\frac{n \pi x}{L}\right)$$

but what if $V$ is not zero? How can I numerically set up the problem to find the eigenvalues and eigenvectors?

We want to solve the eigenproblem $\hat{H}|\psi\rangle=\lambda |\psi\rangle$. Since $|\psi\rangle$ lives in an infinite dimensional Hilbert space, this isn't fit to put on to a computer. Instead, we discretize space into $N$ points. $\vec{\psi}$ is now an $N$ dimensional column vector, and $\bf H$ is now an $N\times N$ matrix. For this post I will write matrices in bold face ($\bf V$), column vectors with an arrow over them ($\vec{\psi}$), and operators with hats over them ($\hat{x}$).

For this post I'll choose to discretize the region from $0$ to $L$ into $N$ points, starting with a point at $x_1=\frac{1}{N+1}L$ and ending at a point $x_N=\frac{N}{N+1}L$. I do this because $\psi(0)=\psi(L)=0$, and I want $N$ nonzero points. Then $x_i=\frac{i}{N+1}L$ and we have our column vector $\psi_i=\psi(x_i)$. Call the step size between these points $b=\frac{L}{N+1}$. There are a few parts to this: How to write $\psi$ as a column vector, how to write the potential $V$ as a matrix, how to write the derivative operator as a matrix, and how to write the full Hamiltonian as a matrix.

Psi as a column vector

Start with an example and take $L=1$. Suppose we want to represent $\psi(x)=\sqrt{\frac{2}{L}}\sin\left(\frac{\pi x}{L}\right)$ as a column vector. If $N=5$, the column vector representing this would be

$$\vec{\psi}=\left[\begin{array}{c} \frac{1}{\sqrt{2}}\\ \sqrt{\frac{3}{2}}\\ \sqrt{2}\\ \sqrt{\frac{3}{2}}\\ \frac{1}{\sqrt{2}} \end{array}\right]$$ This column vector has entries $\psi_i=\psi(x_i)$. Note that the normalization condition $\int_0^L \psi(x)^2dx=1$ turns into the condition $\vec{\psi}^T\vec{\psi} b=1$. (Check it: $\vec{\psi}^T\vec{\psi}=\frac{1}{2}+\frac{3}{2}+2+\frac{3}{2}+\frac{1}{2}=6$ and $b=1/6$. "b" takes the place of "dx"). It's a coincidence that this works exactly for the function we chose, but it tells us we should choose $\vec{\psi}^T\vec{\psi} b=1$ as a normalization condition. If $\psi$ were complex then this should indeed be $\vec{\psi}^\dagger \vec{\psi}b=1$.

V as a matrix

The operator $\hat{V}$ sends $\psi(x)$ to $V(x)\psi(x)$ for some function $V(x)$. So the matrix should send $\psi_i=\psi(x_i)$ to $V(x_i)\psi(x_i)$. This can be represented by a diagonal matrix. For example, with $N=3$,

$${\bf V}=\left[ \begin{array}{ccc} V(x_1)&0&0\\ 0&V(x_2)&0\\ 0&0&V(x_3) \end{array}\right]$$

because

$${\bf V}\vec{\psi}=\left[\begin{array}{c} V(x_1)\psi(x_1)\\ V(x_2)\psi(x_2)\\ V(x_3)\psi(x_3) \end{array}\right]$$

In Mathematica this can be constructed using the Array command (I use a lower case "n" because upper case N is taken):

DiagonalMatrix[Array[V[#/(n + 1) L] &, n]]


The derivative as a matrix

There are many possible ways to approximate the derivative. For example, we could recall the limit definition of the derivative, and send $\psi(x_i)$ to $\frac{1}{b}\left(\psi(x_{i+1})-\psi(x_i)\right)$. The approximation I'll use for the second derivative will be sending $\psi(x_i)$ to $\frac{1}{b^2}\left(\psi(x_{i+1})+\psi(x_{i-1})-2\psi(x_i)\right)$. Let's consider the $N=5$ case again. Then the operator $\frac{\partial^2}{\partial x^2}$ can be written as:

$${\bf D^2}={\bf \frac{\partial^2}{\partial x^2}}=\frac{1}{b^2}\left[ \begin{array}{ccccc} -2 &1&0&0&0\\ 1&-2&1&0&0\\ 0&1&-2&1&0\\ 0&0&1&-2&1\\ 0&0&0&1&-2 \end{array}\right]$$

(I'm giving it the name ${\bf D^2}$, but please don't think it's the square of an actual matrix ${\bf D}$!)

Note that we run into trouble at the corners, but are saved because $\psi(x_0)=\psi(0)=0$, so we can leave terms out. If the boundary conditions were different, we wouldn't be able to represent the derivative as a matrix and would have to add additional terms. So the nice boundary conditions mean we dodged a bullet!

This stripe of diagonal and immediately off-diagonal terms comes up incredibly often in numerics. It's called the "discrete laplacian", and is an example of a tridiagonal matrix. I learned about these from a book A First Course in Computational Physics by DeVries and Hasbun, and I highly recommend it.

The eigenvectors of $\bf D^2$ are the discretized versions of the sine operator, just as how the eigenfunctions of the second derivative operator are sines and cosines.

In Mathematica, you can use the Band[] function to specify these off-diagonal terms:

D2=1/b^2 SparseArray[{Band[{1, 1}] -> -2, Band[{2, 1}] -> 1, Band[{1, 2}] -> 1}, {n,n}]


The total Hamiltonian

We can now write ${\bf H}=-\frac{\hbar^2}{2m} {\bf D^2}+{\bf V}$. For example, for $N=5$, if I put all the Mathematica code together:

$${\bf H}=\left( \begin{array}{ccccc} \frac{\hbar ^2}{b^2 m}+V\left(\frac{L}{6}\right) & -\frac{\hbar ^2}{2 b^2 m} & 0 & 0 & 0 \\ -\frac{\hbar ^2}{2 b^2 m} & \frac{\hbar ^2}{b^2 m}+V\left(\frac{L}{3}\right) & -\frac{\hbar ^2}{2 b^2 m} & 0 & 0 \\ 0 & -\frac{\hbar ^2}{2 b^2 m} & \frac{\hbar ^2}{b^2 m}+V\left(\frac{L}{2}\right) & -\frac{\hbar ^2}{2 b^2 m} & 0 \\ 0 & 0 & -\frac{\hbar ^2}{2 b^2 m} & \frac{\hbar ^2}{b^2 m}+V\left(\frac{2 L}{3}\right) & -\frac{\hbar ^2}{2 b^2 m} \\ 0 & 0 & 0 & -\frac{\hbar ^2}{2 b^2 m} & \frac{\hbar ^2}{b^2 m}+V\left(\frac{5 L}{6}\right) \\ \end{array} \right)$$

n=5;
x[i_]=i L/(n+1);
Vmatrix=DiagonalMatrix[Array[V[x[#]]&,n]];
D2=1/b^2 SparseArray[{Band[{1,1}]->-2,Band[{2,1}]->1,Band[{1,2}]->1},{n,n}];
H=-hbar^2/(2m)D2+Vmatrix


Finding the Eigensystem

It's not easy to find the eigenvalues of a general matrix. Tridiagonal matrices are easier, but not trivial by any means. You will probably want to use a library or consult more advanced references for this! Fortunately, Mathematica has a nice self-contained function for this called Eigensystem[]. The eigenvectors are our orthonormal wavefunctions, and the eigenvalues are our energy levels. Here's a working, full example, with a potential going like $V(x)=1-4(x-L/2)^2$ and 100 points. (Note: this means $H$ is a 100x100 matrix!) Output:

Input:

n=100; L=1; hbar=1; m=1; b=L/(n+1);
V[x_]:=200(1-4(x-L/2)^2);
x[i_]=i L/(n+1);
normalize[v_]:=v/Sqrt[b v.v]; (* ensure that normalize[v].normalize[v]*b==1 *)
Vmatrix=DiagonalMatrix[Array[V[x[#]]&,n]];
D2=1/b^2 SparseArray[{Band[{1,1}]->-2,Band[{2,1}]->1,Band[{1,2}]->1},{n,n}];
H=-hbar^2/(2m)D2+Vmatrix;
eigensystem=Eigensystem[N[H]]; (* find its eigensystem *)
eigensystem={eigensystem[[1]],normalize/@eigensystem[[2]]}; (* normalize all the eigenvectors *)
Table[ListLinePlot[eigensystem[[2]][[-k]],AxesLabel->{"i","\[Psi]"},PlotLabel->"Eigenstate "<>ToString[k]<>", energy="<>ToString[NumberForm[eigensystem[[1]][[-k]],3]]],{k,1,5}]


Setting V=0 and comparing our numerical accuracy, we find the first eight eigenvalues to be:

Energies = {4.9344, 19.7328, 44.381, 78.855, 123.122, 177.138, 240.852, 314.2011}


whereas the exact values are $n^2 \pi^2/2$, giving:

Energies = {4.9348, 19.7392, 44.4132, 78.9568, 123.37, 177.653, 241.805, 315.827}


So this approach is pretty accurate, giving an error of .5% for the eighth energy level.

• Very nice answer, but IMO it might be better to use Pseudocode since not everyone knows what Mathematica syntax means. – Kyle Kanos May 12 '17 at 10:08