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$,

    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:

enter image description here

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

        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:

       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

      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.

  • $\begingroup$ 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. $\endgroup$
    – Alx
    Nov 13, 2019 at 2:12
  • $\begingroup$ edited to answer your questions, L is the first operator I defined and vals was supposed to be hvals. $\endgroup$
    – Misha
    Nov 13, 2019 at 19:30

1 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:

   Evaluate[(Lh[hfuns[[i]], x]) - hvals[[i]]*hfuns[[i]][x]],
    {x, 0, π}], {i, 6}], 3]

enter image description here

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:

    L[Exp[-h[#]/2] hfuns[[i]][#] &, x] - 
    hvals[[i]]*Exp[-h[x]/2] hfuns[[i]][x]], {x, 0, π}], {i, 6}], 3]

and we get:

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.