Worksheet 4, PHY301 - Fall 2003 - Due Nov 14th

Write a C++ code to solve the following eigenvalue problem from quantum mechanics. Use the iterative power method code you developed in the last worksheet.

But first, write and test a new C++ function that computes both the largest and the smallest eigenvalue of a matrix by (1) calling your previous function to find the eigenvalue of largest magnitude; (2) subtracting that eigenvalue times the identity matrix from the original matrix and calling your previous function to find a new extreme eigenvalue; (3) calculating the other extreme eigenvalue of the original matrix by combining these results.

Problem

Find the lowest eigenvalue and eigenvector of a particle in the potential (i.e. solve $(-d^2/dx^2+V(x)) \psi(x) = E \psi(x)$)

V(x) = x² + b x⁴, (1)

for a range of values of b, and make a plot of your result. (For no good reason, this equation is written in units where hbar=2*M=k/2=1. You may well prefer to take hbar=M=1, so that the harmonic oscillator has its familiar ground-state energy of 1/2.)

To do this calculation, you will make the continuous problem into a discrete one by limiting the range in x to -L < x < +L, and representing the wave function by its value at uniformly spaced lattice points within that interval. You will of course have to check that your final results are independent of L and the lattice spacing. (This approach has an enormous range of applications in modern physics and engineering.)

You will need to set up a matrix that represents the operator d²/dx² by finite differences. Use Taylor's series for the wave function expanded about each lattice point to see how to do this. The simplest method uses only the nearest neighbors. However, it is not difficult to figure out how to do better by also using next-nearest neighbors; or also next-next-nearest.

Also compare your result to the result expected expected from first order perturbation theory, which you can calculate with the help of Mathematica. (Recall that for a harmonic oscillator $E=(n+1/2) \hbar \omega$, and the ground-state wavefunction is a Gaussian.)

Also compare your result to the result of a variational calculation using a Gaussian wave function. Explicitly, use a trail wave function of the form psi(x) = c1 * exp(-c2*x^2), where c1 is a normalization constant and c2 is chosen to minimize the expectation value of the Hamiltonian. This calculation gives an upper bound to the ground state energy, and is exact in the limit b->0 since the true wave function for the harmonic oscillator is a Gaussian. Use Mathematica to do this calculation.

(A different approach to converting this Schroedinger Equation problem to matrix form would be to use harmonic oscillator basis states, i.e., the eigenstates of the Hamiltonian with b set to 0. The resulting infinite dimensional matrix would be truncated at some appropriate dimension. Feel free to try it this way also, if you are sufficiently familiar with the raising and lowering operators! (Hint: to find the infinite-dimensional matrix of x^4, just find the infinite-dimensional matrix of x^2 and square it.))