Eigenvalues and numerical eigenfunctions for similar differential operators

Question

I am looking to numerically approximate the eigenvalues and eigenfunctions for a differential operator I am working with, assuming $\pi$ periodic boundary conditions.

Namely, I define the function $h(x) = -ax + \frac{b}{2}\cos(2x)$, with $a$ and $b$ being model parameters, here we can just use $a=1=b$ for now.

Then I consider the differential operator $$Lf(x) = f''(x) + h'(x)f'(x) = e^{-h(x)}\bigg(e^{h(x)}f'(x)\bigg)'$$ with the periodic boundary conditions $f(0) = f(\pi)$. This is the operator whose eigenvalues and eigenfunctions I want. This operator is not self-adjoint with respect to the usual $L^2([0,\pi])$ inner product. However, this operator is self adjoint with respect to the inner product $$(f, g)_h = \int_0^\pi f(x)g(x)e^{h(x)}dx$$ As can be easily shown either numerically or by integration by parts.

In order to get the eigenvalues from NDEigensystem correctly, I define $$T_h: C^\infty(S^1) \to C^\infty(S^1), \quad T_h(f) = e^{h(x)/2}f(x)$$ and notice that $T_{-h} = T_h^{-1}$, and that $(f,g)_h = (T_hf,T_hg)$, where $(f,g)$ is the usual $L^2$ inner product.

I then define

$$L_h: C^\infty(S^1) \to C^\infty(S^1),\quad L_hf := T_hLT_h^{-1}f$$

and observe that $L_h$ is self-adjoint with the unweighted inner product. Thus, if I can get the eigenvalues for $L_h$, I have the eigenvalues for $L$, as eigenvalues for similar linear differential operators are the same.

In order to work with $L_h$, we observe that by wolfram alpha or just long computation, $$\begin{align*} L_hf(x) &= T_hL[e^{−h(x)/2}f(x)] = e^{h(x)/2} e^{-h(x)}\bigg(e^{h(x)}(e^{-h(x)/2}f(x))'\bigg)'\\ &= e^{-h(x)/2}\bigg(e^{h(x)}(e^{-h(x)/2}f(x))'\bigg)'\\ &= h'e^{h(x)/2}(e^{-h(x)/2}f(x))' + e^{h(x)/2}(e^{-h(x)/2}f(x))'' \\ &= f''(x) + f(x)\bigg(-\frac{1}{2}h''(x)-\frac{1}{4}(h'(x))^2 \bigg) \end{align*} $$

And so I use NDeigensystem to get the eigensystem for $L_h$,

    ClearAll
    b:= 1
    c:= 1
    Lh[l_, x_] := D[l[x], x, x] +(-D[h[x],x,x]/2-((D[h[x],x])^2)/4) l[x]
    {hvals, hfuns} = 
    NDEigensystem[{Lh[l, x], l[0] == l[Pi]}, l, {x, 0, Pi}, 6, 
        Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions" \
        -> {"MaxCellMeasure" -> 0.001}}}}]

With the output:

And finally I check that this is indeed the correct eigensystem, via the method described in the answer to my question here: Checking NDEigensystem Results

    Do[Print[Plot[
        Evaluate[Lh[hfuns[[i]], x] - hvals[[i]]*hfuns[[i]][x]], {x, 
        0, \[Pi]}]], {i, 6}]

With the output nicely showing that the difference between $L_hf$ and $\lambda_f$ for $f$ an eigenfunction with corresponding eigenvalue $\lambda_f$ is negligable.

Now, my issue is checking that these are the correct eigenvalues for my original operator $L$ $$Lf(x) = f''(x) + h'(x)f'(x) $$

    L[l_,x_] := D[l[x],x,x] + (D[h[x],x])(D[l[x],x])

As $L_hf = T_h L T_{-h} f$, $f$ an eigenfunction for $L_h$ with eigenvalue $\lambda_f$ implies that $e^{-h(x)/2}f(x)$ is an eigenfunction for $L$ with eigenvalue still $\lambda_f$. However, I am running into issues using mathematica to verify this:

    Do[Print[Plot[
       Evaluate[L[Exp[-h[x]/2]hfuns[[i]], x] -  hvals[[i]]*Exp[-h[x]/2]hfuns[[i]][x]], 
       {x,0, \[Pi]}]], {i, 6}]

returns only axes with no plots, and

    Do[Print[Plot[
      Evaluate[L[FunctionInterpolation[Exp[-h[x]/2]hfuns[[i]],{x,0,Pi}], x] -    hvals[[i]]*FunctionInterpolation[Exp[-h[x]/2]hfuns[[i]],{x,0,Pi}][x]], 
      {x,0, \[Pi]}]], {i, 6}]

Returns errors, although I'm not confident I'm using it correctly.

How should I verify that the eigenvalues and eigenfunctions I'm using are in fact the correct eigensystem for $L$?

Thank you, and sorry for the long post, I wanted to make sure the details were all there.

Answer 1

This is not complete answer, I only show how to overcome your issues with Plot.

Your first check

Do[Print[Plot[ Evaluate[Lh[hfuns[[i]], x] - hvals[[i]]*hfuns[[i]][x]], {x, 0, π}]], {i, 6}]

can be done in the following way:

GraphicsGrid@
 Partition[
  Table[Plot[
   Evaluate[(Lh[hfuns[[i]], x]) - hvals[[i]]*hfuns[[i]][x]],
    {x, 0, π}], {i, 6}], 3]

Now to you operator L:

L[l_, x_] := D[l[x], {x, 2}] + (D[h[x], x]) (D[l[x], x])

To make it work in plotting commands one needs to represent it as pure function:

GraphicsGrid@
 Partition[
  Table[Plot[
   Evaluate[
    L[Exp[-h[#]/2] hfuns[[i]][#] &, x] - 
    hvals[[i]]*Exp[-h[x]/2] hfuns[[i]][x]], {x, 0, π}], {i, 6}], 3]

and we get: