NDEigensystem potential with a singularity

Question

Consider the following eigenvalue equation

$$-\frac{d^2}{dx^2}\Psi_n(x)+\left(gx^4+\frac{1}{x^2}\right)\Psi_n(x)=E_n\Psi_n(x),\qquad x\in(-10,10),\qquad\Psi(-10)=\Psi(10)=0$$

The boundary of $x$ is arbitrary, Ideally it should be $x\in(-\infty,\infty)$ but since we are doing numerics we cant impose that. I want to solve this system using NDEigensystem. I have tried the following so far

n = 25;
l = 10;
g = 1/4;
V[x_] := (g*x^4 + 1/x^2);
op = -D[D[\[Psi][x], x], x] + V[x]*\[Psi][x];
bc = DirichletCondition[\[Psi][x] == 0, True];
{vals, funs} = NDEigensystem[{op, bc}, \[Psi][x], {x, -l, l}, n];

When I run this code, I get all sorts of crazy eigenfunctions with singularities. Not ideal at all. But when I run it with a smaller $l$, say l=Pi, I get much well behaved functions. How should I adjust my boundaries, if ideally they need to be infinity. Should I somehow consider the potential too?

Answer 1

There are some considerations that should be made before solving the equation.

First of all, your potential, $V(x)=gx^4+\frac{1}{x^2}$, is singular at the origin. Hence, an additional condition is needed if you want well-behaved eigenfunctions: vanishing function at $x=0$.

Secondly, $V(x)$ is even: $V(x)=V(-x)$. Therefore, it is sufficient to solve the equation on the semi-infinite interval $[0,\infty)$ with boundary conditions $\Psi_n(0)=\Psi_n(\infty)=0$. However, in practice, this is not possible as you mention correctly. Thus, you have to truncate your domain and solve on $[0,l]$, with $l$ sufficiently large.

The natural question: What does it mean sufficiently large $l$?

In this particular example, asymptotic analysis provides the answer. It can be easily demonstrated that any eigenfunction at large $|x|$ behaves like

$$ \Psi_n(x)\sim \exp\left(-\frac{1}{6}|x|^3\ +\ ...\right)\ ,\qquad (g=1/4)$$

see Bender & Orszag for details on a general setting. It means that eigenfunctions decay extremely fast once $|x|$ is large. For example, $l=6$ leads to $|\Psi_n(l)|\sim 10^{-16}$.

Based on your code, I propose the following: replace {x, -l, l} by {x, 0, l} and set l=10. Thus,

n = 25;
l = 10;
g = 1/4;
V[x_] := (g*x^4 + 1/x^2);
op = -D[D[\[Psi][x], x], x] + V[x]*\[Psi][x];
bc = DirichletCondition[\[Psi][x] == 0, True];
{vals, funs} = NDEigensystem[{op, bc}, \[Psi][x], {x,0, l}, n, Method -> {"SpatialDiscretization" -> {"FiniteElement", 
"MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}}];

I inserted some options in Method to increase the accuracy in calculations. Playing a little bit with the value of $l$, it seems that the first 10 eigenvalues are reliable when $l=10$ and the corresponding eigenfunctions are well-behaved, as expected. See below.