Solution to eigenvalue BVP using NDEigensystem to high precision

Question

I'm trying to solve linear (non-self-adjoint) boundary-value problems to as high precision as possible (optimally 1e-15). For example, the below code solves for the first 5 eigenvalues of the harmonic oscillator,

$$ -u'' + V(x) u = Eu $$

where $V(x) = x^2$ which should be solved subject to $u \to 0$ as $|x| \to \pm\infty$. It then outputs the error compared to the exact values of

$$E_n = (2n + 1), \quad n = 0, 1, \ldots $$

xmax = 100;
numeigs = 5;
maxcell = 0.01;
V[x_] := x^2
Leq = -u''[x] + V[x]*u[x];
{vals, funs} = 
  NDEigensystem[Leq, u[x], {x, -xmax, xmax}, numeigs, 
   Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {MaxCellMeasure -> maxcell}}}}];
vals - Table[(2*n + 1), {n, 0, numeigs - 1}]

with errors:

{2.64357*10^-11, 1.82258*10^-10, 6.51789*10^-10, 1.64085*10^-9, 
 3.36012*10^-9}

There are two key parameters: the xmax condition should tend to infinity and the maxcell should tend to zero. I've tried raising and lowering both but I can't seem to get better than 1e-10 at the first eigenvalue.

I've also tried adding something like this:

arnoldcond := 
  "Eigensystem" -> {"Arnoldi", "MaxIterations" -> Infinity, 
    Tolerance -> 10^(-20)};

but it doesn't seem to help much and to be honest I have no idea what that "Tolerance" actually does.

Additional reference [1] Gaining precision/accuracy with NDEigenvalues

Answer 1

Not a solution using NEigensystem, but I have a package for numerically calculating solutions of eigenvalue problems using the Evans function via the method of compound matrices, which is hosted on github. See my answers to other questions or the github for some more details.

First we install the package (only need to do this the first time):

Needs["PacletManager`"]
PacletInstall["CompoundMatrixMethod", 
    "Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]

Then we first need to turn the ODEs into a matrix form $\mathbf{y}'=\mathbf{A} \cdot \mathbf{y}$, using my function ToMatrixSystem.

Needs["CompoundMatrixMethod`"]
sys = ToMatrixSystem[-u''[x] + V[x] u[x] == e u[x], {u[-L] == 0, 
    u[L] == 0}, {u}, {x, -L, L}, e] /. V -> Function[{x}, x^2]

Now the function Evans will calculate the Evans function (also known as the Miss-Distance function) for any given value of $\lambda$; this is an analytic function whose roots coincide with eigenvalues of the original equation.

Plot[Evans[e, sys /. L -> 5], {e, 0, 10}]

Now FindRoot will find roots, and you can up the precision as required:

(3 - e) /. 
  FindRoot[Evans[e, sys /. L -> 10, NormalizationConstants -> 0, 
    WorkingPrecision -> 50], {e, 2}, WorkingPrecision -> 50] // Quiet
(* -1.477961860221512919729879633184*10^-19 *)

As this system is even, you can find even modes with only one infinite endpoint:

sys2 = ToMatrixSystem[-u''[x] + V[x] u[x] == e u[x], {u'[0] == 0, 
    u[L] == 0}, {u}, {x, 0, L}, e] /. V -> Function[{x}, x^2]

Now the default normalisation included in Evans works better, but it misses the odd eigenvalues:

Plot[Evans[e, sys2 /. L -> 5], {e, 0, 20}]

(5 - e) /. 
  FindRoot[Evans[e, sys2 /. L -> 10, WorkingPrecision -> 50], {e, 
    4.01}, WorkingPrecision -> 50] // Quiet
(* 4.73068503328772292182316234385043*10^-17 *)

Answer 2

The finite element code that NDEigensystem uses is available only in machine precision (in Version 12.0). So in this case, if mesh refinement does not bring further improvement, you are out of luck with the FEM.