Solving the elliptic eigenproblem with periodic boundary conditions

Question

Given an energy level $\mu$, I'm looking to calculate the eigenvectors corresponding to the time-independent Schrodinger operator on the torus (that is, periodic boundary conditions)-- $H = -h^2 \Delta + V$ for some nice, bounded potential $V$ -- whose eigenvalues are nearest to $\mu$.

I was wondering if there are some standard techniques for solving this particular problem, where $h$ can be rather small, say $h = 0.001$.

My current method is the naive approach (by my reckoning): I use the DFT to convert into phase space, whereby I can compute the map $\hat{u} \mapsto \hat{(H u)} = h^2 |k|^2 \hat{u} + \hat{V} \ast \hat{u}$, and then proceed by using the standard Krylov method for the matrix $(\hat{H} - \mu I)^{-1}$, where the map corresponding to that inverse is computed using biconjugate gradients.

There are several key problems with this approach. The first is that $h$ is very small, so for certain energy levels I must use the DFT up to a very high frequency in order to get an accurate discretization, but this means that $\hat{H}$ is a both large (larger than $10000\times 10000$) and dense, so Krylov methods are still quite expensive, especially since on each Arnoldi iteration I have to do a full set CG iterations.

Are there any other standard methods that I should try? I'd be grateful for any suggestions. Even if it's something obvious, I probably haven't thought of it yet because I'm not very well-versed in numerical eigenvalue problems, especially not in the case of differential operators.

Answer 1

A question first: What program are you using to perform the eigendecomposition of your matrix?

Rather than look for an answer in phase space, you may be better off considering your problem in physical space instead. In that case, you can use the standard finite difference formula to approximate $\Delta u$ at each point, then form a matrix corresponding to that system, etc.

The reason for this is that the resulting matrix is sparse, so there is less arithmetic to do to compute a matrix-vector product. If you're looking at high frequencies, you'll need a lot of points to resolve them and thus your system will be of high dimension; sparsity really saves you here. Matlab and Scipy both have built-in sparse matrix data types, which you'll have to use in order to take advantage of this.

You also mentioned using the Arnoldi process. Using the Lanczos process will be faster, since your problem is symmetric. If you're using a program smart enough to detect that your matrix is symmetric and act accordingly then you needn't worry about this, but if you explicitly selected Arnoldi then you could be doing unnecessary work solving upper Hessenberg systems when it could be solving tridiagonal systems. Nonetheless, you're still finding the eigenvalues closest to zero of the operator $(H-\mu I)^{-1}$ -- that at least doesn't change.

I would also point you to this thread as a caveat: eigensolvers can give erroneous results in a way which, if you're interpreting the results as describing a quantum system, are very unphysical. Long story short, don't look for the whole spectrum of a matrix, just the first few eigenvalues. It sounds like that's what you're doing anyway.