Checking NDEigensystem Results

Question

I'm looking to verify the output of a call to NDEigensystem. I'm doing this by plotting the operator acting on the Interpolating Function outputs versus the eigenvalue times the interpolating eigenfunctions. Something is going wrong, as it seems to me that by the properties of eigenvalues and eigenvectors, the plots should line up.

Here's my code:

ClearAll
b:= 1
c:= 1

h[x_] := -b x + c Cos[2x] / 2

Lh[l_, x_] := (D[l[z],z,z] -(1/4)(D[h[z],z])^2l[z]) /. z :> x

{vals, funs} = NDEigensystem[{-Lh[l,x], l[0] == l[Pi]}, l[x], {x, 0, Pi}, 5]

For[i=1, i<6, i++, f[x_] := funs[[i]]; v := vals[[i]]; 
P1 =  Plot[-Lh[f,x], {x,0, Pi}, PlotStyle -> Red ]; 
P2 = Plot[v f[x], {x,0,Pi}, Axes->False];
Print[Overlay[{P1, P2}]]]

And my ouput is:

For $L$ a differential operator on $L^2(S^1)$ with eigenfunction $u$ and eigenvector $\lambda$, we should have $Lu(x) = \lambda u(x)\, \forall x \in S^1$. However, according to the plots, this seems to not be the case.

Where is my error? Thanks!

Answer 1

I think something is not quite right with way you process the equations. If you use this you get what you expect:

b = 1;
c = 1;
h[x_] = -b x + c Cos[2 x]/2;
Lh[l_, x_] := D[l[x], x, x] - (1/4) (D[h[x], x])^2 l[x]
{vals, funs} = 
  NDEigensystem[{-Lh[l, x], l[0] == l[Pi]}, l, {x, 0, Pi}, 6, 
   Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions" \
-> {"MaxCellMeasure" -> 0.001}}}}];
Do[Print[Plot[
   Evaluate[-Lh[funs[[i]], x] - vals[[i]]*funs[[i]][x]], {x, 
    0, \[Pi]}]], {i, 6}]

If you need a better resolution, then you'd need to refine the mesh. One other point that is important: With your approach, you are also testing the ability to find higher order derivatives of interpolating functions. That in it self is prone to numerical issues. So you are testing two things.

Answer 2

For NDEigensystem[], only homogeneous boundary conditions can be set. In this case, we have

b = 1;
c = 1;

h[x_] := -b x + c Cos[2 x]/2

eq = -(D[l[z], z, z] - (1/4) (D[h[z], z])^2 l[z]);

{vals, funs} = 
 NDEigensystem[{eq, DirichletCondition[l[z] == 0, z == 0 || z == Pi]},
   l[z], {z, 0, Pi}, 5]

Table[Plot[{eq /. l[z] -> funs[[i]], vals[[i]] funs[[i]]}, {z, 0, 
   Pi}], {i, 1, 5}]