# Eigenvalues and numerical eigenfunctions for similar differential operators

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.

• Your Do[Print[Plot[Evaluate[Lh[hfuns[[i]], x] - hvals[[i]]*hfuns[[i]][x]], {x, 0, \[Pi]}]], {i, 6}] can be presented as GraphicsGrid@Partition[Table[Plot[Evaluate[(Lh[hfuns[[i]], \[FormalX]] /. \[FormalX] -> x) - hvals[[i]]*hfuns[[i]][x]], {x, 0, \[Pi]}], {i, 6}], 3]. And for your two last commands: what are L, vals? They are not defined previously in your question. – Alx Nov 13 '19 at 2:12
• edited to answer your questions, L is the first operator I defined and vals was supposed to be hvals. – Misha Nov 13 '19 at 19:30

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

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: