NDEigenvalues complains not Hermitian with large dimension differential operator

Question

The following snippet calculates eigenvalues and eigenfunctions of a null operator (just as an example):

n = 10;
us = Table[Symbol[StringJoin["u", ToString[i]]][x, y], {i, n}];
NDEigenvalues[
 Join[
  Table[0.0, {i, n}],
  Table[DirichletCondition[us[[i]] == 0, True], {i, n}]
  ],
 us,
 {x, y} \[Element] Rectangle[{-5.0, -5.0}, {5.0, 5.0}],
 6
 ]

With n less than 10, everything works perfectly. For n>10, NDEigenvalues complains Eigensystem::herm: The matrix SparseArray[Automatic,<<3>>] is not Hermitian or real and symmetric.

Is this a precision issue? How to avoid this?

Edit: The official suggestion to avoid this problem is to set InterpolationOrder, as shown in the link posted by @user21.

Answer 1

Here is how you can avoid that:

n = 10;
us = Table[Symbol[StringJoin["u",
    StringPadLeft[ToString[i], 2, "0"]]][t, x, y], {i, n}];
NDEigenvalues[
 Join[Table[0.0, {i, n}], 
  Table[DirichletCondition[us[[i]] == 0, True], {i, n}]], us, {x, 
   y} \[Element] Rectangle[{-5.0, -5.0}, {5.0, 5.0}], 6]
{1.5777218104420236`*^-29, 1.8932661725304283`*^-29, \
2.208810534618833`*^-29, 2.8398992587956425`*^-29, \
3.470987982972452`*^-29, 5.048709793414476`*^-29}

Note that I changed the way the dependent variables are generated to:

us = Table[Symbol[StringJoin["u",
StringPadLeft[ToString[i], 2, "0"]]][t, x, y], {i, n}]
{u01[x, y], u02[x, y], u03[x, y], u04[x, y], u05[x, y], u06[x, y], 
 u07[x, y], u08[x, y], u09[x, y], u10[x, y]}

This has the property that

Sort[us]==us
True

In essence what you are seeing is an issue discussed section Ordering of Dependent Variable Names in the documentation; though it might not be quite that obvious.

How did I arrive at this conclusion? One thing you can do is to follow this post to see what is happening inside. We start by making a coarser mesh and using ProcessEquations to generate the FEM data.

Needs["NDSolve`FEM`"]
m = ToElementMesh[Rectangle[{-5.0, -5.0}, {5.0, 5.0}], 
   "MaxCellMeasure" -> 20, "MeshOrder" -> 1];
n = 10;
us = Table[Symbol[StringJoin["u", ToString[i]]][t, x, y], {i, n}];
(*us = Table[Symbol[StringJoin["u",
    StringPadLeft[ToString[i], 2, "0"]]][t, x, y], {i, n}];*)
{state} = 
  NDSolve`ProcessEquations[
   Join[Table[D[us[[i]], t] == 0, {i, n}], 
    Table[us[[i]] == 0 /. t -> 0, {i, n}]], 
   us, {t, 0, 1}, {x, y} \[Element] m, 
   Method -> {"MethodOfLines", "TemporalVariable" -> t, 
     "SpatialDiscretization" -> {"FiniteElement"}}];

Next we extract the FEM data from the NDSolve state object.

femd = state["FiniteElementData"];
bcd = femd["BoundaryConditionData"];
md = femd["FEMMethodData"];
pdec = femd["PDECoefficientData"];
vd = state["VariableData"];
sd = state["SolutionData"][[1]];

And discretize the system:

{l, s, d, mm} = DiscretizePDE[pdec, md, sd]["SystemMatrices"];

Now, look at that:

MatrixPlot[d]

This does not look, right. It seems that the ordering of the variables is not correct. You can try with n=9 to see the difference. This gave me clue to look at the dependent variable ordering. And sure enough looking at

NDSolve`SolutionDataComponent[vd, "DependentVariables"]
{u1, u10, u2, u3, u4, u5, u6, u7, u8, u9}

revealed the culprit. And the fix is, as described above, to use a dependent variable ordering that does not suffer from an internal reordering. Surely this is very unfortunate behavior of the NDSolve function family in conjunction to the FEM. The fact that NDSolve did this reordering was done long before the FEM got implemented and it seems very hard to fix. Which is unfortunate. Sorry about that.