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, 

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


2 Answers 2


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):

    "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.

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}]

enter image description here

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}]

enter image description here

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

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.