Activity: Kinetic energy

Computational Physics Lab II 2022
Students implement a finite-difference approximation for the kinetic energy operator as a matrix, and then use numpy to solve for eigenvalues and eigenstates, which they visualize.
What students learn
  1. The wave function relates to a ket by \(\psi(x) = \langle x|\psi\rangle\).
  2. A wave function can be represented as an array holding the value of the wave function on a grid of points.
  3. In this “position basis," the potential energy operator is diagonal, but the kinetic energy operator is not.
  4. A differential operator like \(\hat T\) can be approximated as a finite difference.

This activity introduces finite difference approximations, follows Position operator.

We have already seens several representations of the state of a particle moving in one dimension. We looked at wave functions and representation in a set of sinusoidal basis functions. The sinusoidal basis set had the advantage that we could use a finite amount of information to describe a smooth function at least approximately.

You might think a wave function can be described with a finite amount of information, but that is only true for a function that can be described analytically. In the general case, we have no such luck, and require an infinite amount of information, essentially the value of the function at every possible position \(x\).

There is another very effective finite representation for wave functions, which you've actually been using for weeks now (at least in effect), but we haven't talked about as a representation. That is the value of the function on a grid of regularly spaced points separated by a distance \(\Delta x\). \begin{align} |\psi\rangle &\,\dot= \begin{pmatrix} \psi(\Delta x) \\ \psi(2\Delta x) \\ \psi(3\Delta x) \\ \psi(4\Delta x) \\ \vdots \\ \psi(L-2\Delta x) \\ \psi(L-\Delta x) \end{pmatrix} \end{align} This “discretized wave function” representation is what you have already been using when you create a plot, and also for numerically integrating to find inner products.

Note: I am not including in this vector \(\psi(0)\) (or \(\psi(L)\) at the other end). This is because the boundary condition at the edge of the box requires that the wave function have a zero value at those two points. We could instead have chosen to include those two points, and then manually forced their values to be zero.

The kinetic energy

The kinetic energy of a particle in one dimension is given by \begin{align} \hat T &= \frac{\hat p^2}{2m} = -\frac{\hbar^2}{2m\Delta x^2}\frac{d^2}{dx^2} \end{align} Now this is a new beast for you, which we call a differential operator. Whenever you are confused by an operator, it helps to operate it on something. \begin{align} \hat T\psi(x) &= -\frac{\hbar^2}{2m}\left.\frac{d^2\psi}{dx^2}\right|_x \end{align} Our question for today is how we can represent this operator in our new “discretized wave function” representation. This requires us to think about what a derivative means, and I'll start with a first derivative: \begin{align} \left.\frac{d\psi}{dx}\right|_x &= \lim_{\Delta x\rightarrow 0}\frac{\psi(x+\Delta x/2) - \psi(x-\Delta x/2)}{\Delta x} \end{align} where I have used a centered difference, because it is symmetric. Now if \(\Delta x\) is reasonably small, we can just omit the limit, which is called a finite difference approximation. Now we want a second derivative, so we need to repeat this. \begin{align} \frac{d^2\psi}{dx^2}(x) &= \lim_{\Delta x\rightarrow 0}\frac{\left.\frac{d\psi}{dx}\right|_{x+\Delta x/2} - \left.\frac{d\psi}{dx}\right|_{x-\Delta x/2}}{\Delta x} \\ &\approx \frac{\left(\psi(x+\Delta x) - \psi(x)\right)-\left(\psi(x)-\psi(x-\Delta x)\right)}{\Delta x^2} \\ &= \frac{\psi(x+\Delta x) +\psi(x-\Delta x)- 2\psi(x)}{\Delta x^2} \end{align} So you can see that the second derivative at \(x\) is sort of a difference between the average of \(\psi\) at the surrounding points and the value of \(\psi(x)\). We can now plug this approximation into our definition for the kinetic energy operator to find: \begin{align} \hat T\psi(x) &\approx \frac{\hbar^2}{2m\Delta x^2}\left(2\psi(x)-\psi(x+\Delta x)-\psi(x-\Delta x)\right) \label{eq:finite-diff} \end{align} This equation is sufficient to express the kinetic energy operator as a matrix in terms of our discretized wave function representation. I will give you the result here, and ask you to prove it in a moment: \begin{align} \hat T &\,\dot= \frac{\hbar^2}{2m\Delta x^2} \begin{pmatrix} 2 & -1 & 0 & 0 & \cdots \\ -1 & 2 & -1 & 0 & \cdots \\ 0 & -1 & 2 & -1 & \cdots \\ 0 & 0 & -1 & 2 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{pmatrix} \label{eq:T-matrix} \\ \begin{pmatrix} \hat T\psi(\Delta x) \\ \hat T\psi(2\Delta x) \\ \hat T\psi(3\Delta x) \\ \hat T\psi(4\Delta x) \\ \vdots \end{pmatrix} &= \frac{\hbar^2}{2m\Delta x^2} \begin{pmatrix} 2 & -1 & 0 & 0 & \cdots \\ -1 & 2 & -1 & 0 & \cdots \\ 0 & -1 & 2 & -1 & \cdots \\ 0 & 0 & -1 & 2 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{pmatrix} \begin{pmatrix} \psi(\Delta x) \\ \psi(2\Delta x) \\ \psi(3\Delta x) \\ \psi(4\Delta x) \\ \vdots \end{pmatrix} \end{align}


  1. On paper (or whiteboard) confirm by matrix multiplication that the matrix in Eq. \ref{eq:T-matrix} is equivalent to the finite difference equation in Eq. \ref{eq:finite-diff}. Raise your hand.
  2. Create a 2D array (or matrix) that represents the kinetic energy operator in a discrete "wave function" representation (i.e. Eq. \ref{eq:T-matrix}).
  3. Find the eigenvalues and eigenfunctions/eigenvectors of the kinetic energy operator. (use numpy)
  4. Plot the 4 eigenfunctions with the lowest energy. Raise your hand.
  5. Create a second plot in which you visualize the eigenvalues.
  6. Make your \(\Delta x\) smaller, and see how these things change.
Extra fun

Construct a matrix for a harmonic potential energy operator \begin{align} \hat V &= \frac12 k \hat x^2 \end{align} This will require you to make use of the fact that \begin{align} \hat x \left|{\psi}\right\rangle \dot=\, x\psi(x) \end{align} which means that \begin{align} \begin{pmatrix} (\Delta x)\psi(\Delta x) \\ (2\Delta x)\psi(2\Delta x) \\ (3\Delta x)\psi(3\Delta x) \\ (4\Delta x)\psi(4\Delta x) \\ \vdots \end{pmatrix} &= \begin{pmatrix} x_{11} & x_{12} & x_{13} & \cdots \\ x_{21} & x_{22} & x_{23} & \cdots \\ x_{31} & x_{32} & x_{33} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix} \begin{pmatrix} \psi(\Delta x) \\ \psi(2\Delta x) \\ \psi(3\Delta x) \\ \psi(4\Delta x) \\ \vdots \end{pmatrix} \end{align} where \(x_{11}\) etc are unknowns that you must determine.

Once you have created a matrix representation for \(\hat x\) and \(\hat V\) in the discretized wavefunction representation, solve for the eigenvectors and eigenvalues of \(\hat V\) and visualize its eigenvectors.

Then on a separate figure plot the potential \(V(x)\) and visualize the eigenvalues as horizontal lines.

Springy fun

Construct a matrix for a Hamiltonian \begin{align} \hat H &= \hat T + \hat V \end{align} Solve for the eigenvectors and eigenvalues, and visualize a few of the lowest energy eigenvalues and eigenvectors.

In a separate figure plot the potential \(V(x)\) and visualize the energy eigenvalues as horizontal lines.

How are your results affected by changing the spring constant \(k\)?

Learning Outcomes
  • ph366: 1) Write functions and entire programs in python
  • ph366: 2) Apply the python programming language to solve scientific problems
  • ph366: 3) Use the matplotlib and numpy packages
  • ph366: 4) Model the physical systems studied in the course