Computing the eigenvectors $\psi(x)$ and eigenvalues $E$ of some Hamiltonian $H$ belongs to the key problems of quantum mechanics. For the one-dimensional Schrödinger equation we have to solve

\begin{equation}

E\psi(x) = H\psi(x) = -\frac{1}{2}\frac{\mathrm{d}^2\psi(x)}{\mathrm{d} x^2} + V(x)\psi(x)

\end{equation}for some given potential $V(x)$. Analytical solutions, however, are scare goods. Thus, one has to rely on approximations or numerical methods.

In principle, numerical approaches are straight forward. Just replace the Hermitian Hamiltonian by some Hermitian matrix by discretizing the differential operators, e.g., by finite differences or pseudospectral methods. Then utilize the computational routines of your favorite computational linear algebra package to solve the eigenvalue problem. Often these packages are equipped with special routines that can take advantage of hermiticity. Below I give an example, which calculates the eigenvalues and the eigenvectors of the harmonic oscillator with $V(x)=x^2/2$ in Python. Differentials are approximated via 2nd order finite differences. Numerical results for the discrete Hamiltonian matrix are compared to the known eigenvalues and eigenvectors. One clearly sees that the small eigenvalues are well reproduced. The discrete Hamiltonian features, however, also large eigenvalues and eigenstates which are very different from the eigenvalues and eigenvectors of the related contiguous problem.

Our naive approach works reasonably well. It has, however, a serious drawback. It becomes computationally very demanding as soon as one tries to tackle two- or three-dimensional problems because the matrices become easily very large. So large that one can not compute its eigenvalues and eigenvectors in a reasonable amount of time. Even worse, the matrices may become so large that they no longer fit in the computers main memory. Fortunately, one is usually not interested in all eigenvalues, knowing some, e.g., a few of the smallest ones, might be sufficient. In this case the Lanczos algorithm is often the method of choice. This algorithm basically constructs from the original large Hermitian matrix a smaller tridiagonal symmetric matrix that has eigenvalues that are close to some eigenvalues of the larger original matrix. A further nice feature of the Lanczos algorithm is that it does not require to store the matrix explicitly in memory, a routine that computes matrix-vector products of the matrix is sufficient.

#!/usr/bin/env python # import useful modules import matplotlib from math import factorial from numpy import * from pylab import * from numpy.polynomial.hermite import * # use LaTeX, choose nice some looking fonts and tweak some settings matplotlib.rc('font', family='serif') matplotlib.rc('font', size=16) matplotlib.rc('legend', fontsize=16) matplotlib.rc('legend', numpoints=1) matplotlib.rc('legend', handlelength=1.5) matplotlib.rc('legend', frameon=False) matplotlib.rc('xtick.major', pad=7) matplotlib.rc('xtick.minor', pad=7) matplotlib.rc('lines', lw=1.5) matplotlib.rc('text', usetex=True) matplotlib.rc('text.latex', preamble=[r'\usepackage[T1]{fontenc}', r'\usepackage{amsmath}', r'\usepackage{txfonts}', r'\usepackage{textcomp}']) def V(x): return 0.5*x**2 close('all') figure(figsize=(6, 4.5)) # define grid N=512 # number of grid points x0, x1=-10., 10. # grid boundaries dx=(x1-x0)/(N-1) # grid spacing x=linspace(x0, x1, N) # grid points # setup Hamiltonian H_kin=zeros([N, N]) for k in range(0, N): if k>0: H_kin[k-1, k]=-0.5/dx**2 H_kin[k, k]=1./dx**2 if k+1<N: H_kin[k+1, k]=-0.5/dx**2 H_pot=diag(V(x)) H=H_kin+H_pot # compute eigenvalues E=eigvalsh(H) # plot first eigenvalues n=100 plot(array(range(0, n)), array(range(0, n))+0.5, 'rx', mew=1.25, ms=8, label='harmonic oscillator eigenvalues') plot(array(range(0, n)), sort(E)[0:n], '.', mew=1.25, ms=8, mfc='none', mec='b', label='eigenvalues of the finite differences Hamiltonian') xlabel('index') ylabel('energy') legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, borderaxespad=0.) gca().set_position((0.15, 0.15, 0.8, 0.625)) show() savefig('eigenvalues.png') # compute eigenvectors E, psi_E=eigh(H) # plot some eigenvectors clf() for n in [0, 1, 2, 20]: c=zeros(n+1) c[n]=1 plot(x, psi_E[:, n]/sqrt(dx), label='$n=%i$' % n) legend(ncol=2) xlabel(r'$x$') ylabel(r'$\Psi_n(x)$') gca().set_ylim(-0.8, 1.2) gca().set_position((0.15, 0.15, 0.8, 0.75)) show() savefig('eigenvectors.png') # plot error of some eigenvectors clf() for n in [0, 1, 2, 20]: c=zeros(n+1) c[n]=1 plot(x, (psi_E[:, n]/sqrt(dx))**2 - (hermval(x, c)*exp(-0.5*x**2)/sqrt(sqrt(pi)*2**n*factorial(n)))**2, label='$n=%i$' % n) legend(ncol=2) xlabel(r'$x$') ylabel(r'$|\Psi_n(x)|^2-|\Psi_{n,\mathrm{exact}}(x)|^2$') gca().set_ylim(-0.004, 0.006) gca().set_position((0.2, 0.15, 0.75, 0.75)) show() savefig('eigenvectors_error.png')