How do you improve the accuracy of a finite difference method for finding the eigensystem of a singular linear ODE

Question

I am attempting to solve an equation of the type:

$ \left( -\tfrac{\partial^2}{\partial x^2} - f\left(x\right) \right) \psi(x) = \lambda \psi(x) $

Where $f(x)$ has a simple pole at $0$, for the smallest $N$ eigenvalues and eigenvectors. The boundary conditions are: $\psi(0) = 0$ and $\psi(R)=0$, and I'm only looking at the function over $(0,R]$.

However, if I do a very simple, evenly spaced finite difference method, The smallest eigenvalue is very inaccurate, (sometimes there is a "false" eigenvalue that is several orders of magnitude more negative than the one I know should be there, the real "first eigenvalue" becomes the second, but is still poor).

What affects the accuracy of such a finite difference scheme? I assume that the singularity is what is causing the problem, and that an unevenly spaced grid would improve things significantly, are there any papers that can point me towards a good non-uniform finite difference method? But perhaps a higher order difference scheme would improve it more? How do you decide (or is it just "try both and see")

note: my finite difference scheme is symmetric tridiagonal where the 3 diagonals are:

$\left( -\frac{1}{2 \Delta^2}, \; \frac{1}{\Delta^2} - f(x), \; -\frac{1}{2 \Delta^2} \right)$

Where $\Delta$ is the grid spacing. And I am solving the matrix using a direct symmetric solver (I am assuming that the accuracy is not affected drastically by the solver, am I wrong?)

Answer 1

If you want to increase the accuracy of a finite difference scheme, you can always try increasing the degree of your stencil. On equidistant points, though, this can lead to numerical instabilities. To avoid these problems and still get high accuracy, I would suggest using Spectral Methods.

If your problem has fixed poles, you can try to get around them by splitting your domain and solving two coupled problems.

The Chebfun system (disclaimer: of which I am one of the developers), uses the above mentioned techniques and you could give your problem a quick spin in the chebgui interface. I'd try it myself, but I don't know what your domain or $f(x)$ is.

I've attached a screenshot of the chebgui window for computing the first six eigenvalues and eigenmodes of $-u''(x) - u(x)x = \lambda u$ on $[-1,1]$.

Update

If you want to solve this problem without getting too much into Chebfun, all the details should be in Chapter 9 of Nick Trefethen's book "Spectral Methods in Matlab".

Answer 2

One way to quickly make things better (though likely not much better) is to consider the similarity between the lowest order finite difference methods you use and the lowest order finite element method. If you compute the tri-diagonal matrix you get from using linear finite element shape functions in 1d, then the discretization of the second derivatives will look exactly the same (up to a factor $\Delta x$ but you'll get a different term for what comes off $f(x)\psi(x)$. I don't know how $f(x)$ looks in your case, but where right now you use $f(x_i)$, it will instead be something like $\int_{x_{i-1}}^{x_{i+1}} f(x) \varphi_i(x)$ where $\varphi_i(x)$ is the hat function that peaks at $x_i$. If $f(x)$ is sufficiently simple, then you can compute this integral exactly, and it will provide a more accurate matrix of which you have to find the eigenvalue.

Of course, if you already do finite elements, you might as well invest in using higher order elements which aren't that much more difficult in 1d either.