Gaining precision/accuracy with NDEigenvalues

Question

See further down for an important note

Background

I study (one component of) the semi-classical Pauli operator, $$ P_h=-h^2\Delta+ih(-y,x)\cdot\nabla+\frac{x^2+y^2}{4}-h. $$ For this particular question, let us consider the Dirichlet problem in the disk $x^2+y^2<1$ (I work in more general settings, but my question applies to this case). I am interested in the size of the lowest eigenvalue $\lambda_1(h)$ of the Pauli operator (it will be exponentially small, like $\exp(-C/h)$ for some $C>0$), as the semiclassical parameter $h\to 0$. For this purpose I am trying to study $\lambda_1(h)$ numerically for small $h$.

Current code

My code looks at the moment like this:

paulieigenvalue[h_, cells_] := NDEigenvalues[
 {(1/4)*((-4*h + x^2 + y^2)*u[x, y] - 4*h*(I*x*D[u[x, y],y] - 
        I*y*D[u[x, y],x] + h*(D[u[x, y],{x,2}] + D[u[x, y],{y,2}]))), 
    DirichletCondition[u[x, y] == 0, True]}, 
u[x, y], 
Element[{x, y}, Disk[]],
1, 
Method -> {"SpatialDiscretization" -> {
  "FiniteElement", {"MeshOptions" -> {MaxCellMeasure -> cells}}}}
]

The input

paulieigenvalue[0.001, 0.001]

gives $$\bigl\{1.542\times 10^{-6}-3.09081\times 10^{-20}i\bigr\},$$ while if I change MaxCellMeasure to $0.0001$,

paulieigenvalue[0.001, 0.0001]

I get $$\bigl\{1.91835\times 10^{-8} - 7.75138\times 10^{-21} i\bigr\}.$$

Even though I know that the eigenvalues should be real, the tiny imaginary parts is not the big issue for me. The big issue is that the real part changes so much, almost by a factor of $100$, when I change the MaxCellMeasure option. For this reason I cannot trust (any of) these values.

Questions

1) Is there some way to gain precision with NDEigenvalues, (except playing with MaxCellMeasure, or playing with it in a more controlled way)?

2) Is there some way, with the built-in functions, to get some error control?

3) Is there a better way in Mathematica to obtain, with high numerical accuracy, the eigenvalue $\lambda_1(h)$ with high precision for small $h$?

Update: An important note

I realize this is not fair to Mathematica. In fact, by calculating the energy of $$ u(x,y)=\exp((1-x^2+y^2)/(4h))-\exp(-(1-x^2+y^2)/(4h)) $$ one finds the true inequality $$ \lambda_1(h)<\frac{h \left(e^{\left.-\frac{1}{4}\right/h}-4 h \sinh \left(\frac{1}{4 h}\right)\right)}{1-2 h \sinh \left(\frac{1}{2 h}\right)}. $$ Inserting $h=0.0001$ into the right-hand side, one gets the tiny number $$ 3.369\times 10^{-2172}. $$ I certainly need to do some transformation before I can attack, numerically, this problem. (In this case that transformation has typically to scale the operator by a factor of $e^{-1/(2h)}$, but in general I'm a bit lost.) In any case, I let the question stand, since I am still curious about if it is possible to work on the precision/accuracy in the eigenvalue solver of Mathematica (perhaps in less extreme situations).

Answer 1

This seems to be a hard problem to solve. Here are a few thoughs: First, I'd like to point out that NDEigensystem is based on a Finite Element discretization, similar to what was shown here and here. This means that the solution to the eigen problem will depend on the underlying mesh that is a discrete representation to the region specified. So the first thing to try is to see how well the mesh actually approximates the region. An easy way to do that is the following:

Needs["NDSolve`FEM`"]
m = ToElementMesh[Disk[], MaxCellMeasure -> 0.001];
\[Pi] - Total[Flatten[m["MeshElementMeasure"]]]
2.1839725583561176`*^-8

This computes the measure of each element and we then sum those up; for the Disk[] region this should sum up to $\pi$ - it is essentially the same as integrating 1 over the region. More on this in the ElementMesh generation tutorial in the Region Approximation Quality Section.

This gives us a 10^-8 difference in the exact region and the discretized region. So that is the ball park in which one can expect the solution accuracy to be - at best.

Now, I have rewritten your function to take a mesh as input:

paulieigenvalue[h_, m_] := 
 NDEigenvalues[{(1/4)*((-4*h + x^2 + y^2)*u[x, y] - 
      4*h*(I*x*D[u[x, y], y] - I*y*D[u[x, y], x] + 
         h*(D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}]))), 
   DirichletCondition[u[x, y] == 0, True]}, u[x, y], 
  Element[{x, y}, m], 1]

And using this with the just generated mesh we get:

paulieigenvalue[0.001, m]

NDEigenvalues::femcsdp: -- Message text not found -- (4.53756060549885`) (9.11215726082907`) (82.6939316356966`) (2)

{1.63826*10^-6 + 6.06951*10^-20 I}

I am sorry about the missing message text; I just fixed that in the code and the next release will have this fixed. The message should be the same as it's NDSolve cousin:

Message[NDSolve::femcsdp, 4.53756060549885` , 9.11215726082907` , \
82.6939316356966`, 2]

NDSolve::femcsdp: The product of 2 times the Damkoehler number (4.53756060549885`) and the Peclet number (9.11215726082907`) is 82.6939316356966` and is larger than the mesh order (2), and the computed result may not be stable. Adding artificial diffusion may help.

This tells you that this is a convection dominated PDE and there is a tutorial section about Stabilization of Convection-Dominated Equations. These type of problems are hard to solve for the FEM. One is, to some degree, by refining the mesh.

m = ToElementMesh[Disk[], MaxCellMeasure -> 0.0000125]
ElementMesh[{{-1.,1.}, {-1., 1.}}, {TriangleElement["<" 396918 ">"]}]

This mesh has about 400T elements. And this is the amount needed to at least satisfy the mesh requirement but that still does not grantee a good solution.

(* this will requite a bit of memory and a few minutes *)
paulieigenvalue[0.001, m]
{2.9862911638824127`*^-10 - 9.034619159308821`*^-20 I}
Pi - Total[Flatten[m["MeshElementMeasure"]]]
3.5340619319867983`*^-12

You can increase the accuracy of the mesh a bit (and reduce the number of elements - to keep computational cost reasonable)

m = ToElementMesh[Disk[], AccuracyGoal -> 8, 
  MaxCellMeasure -> 0.00002]

Pi - Total[Flatten[m["MeshElementMeasure"]]]
-5.3290705182007514`*^-14

(* this will requite a bit of memory *)
paulieigenvalue[0.001, m]
{7.646825379547697`*^-10 + 3.960491626284256`*^-20 I}

To your questions:

1) Yes, use PrecisionGoal and AccuracyGoal - this may not always work from an NDSolve/NDEigensystem level (this is subject to a future improvement) but it works on the ToElementMesh level.

2) I tried so show that with the post

3) Not easily; as the PDE itself poses trouble and the FEM is implemented for machine precision and uses second order mesh approximations.

A note on why you did not get a warning when you specified the Disk region directly. Ideally, (in a future version) NDEigensystem would generate a mesh that is appropriate for the PDE but in this case that would have resulted in a large mesh any ways.

A further thing to note is that h is an average mesh cell diameter.

I can not exactly predict if and how scaling will affect this - but if you do try that out, I'd be interested in hearing the result. But do think it's worthwhile to try.