Plot eigenvalues of fractional Laplacian in 1D

Question

Let $D=(-1,1)$ and $\alpha \in(0,2)$.

Question: How can I plot the eigenvalues $(\lambda_n)_{n \in \mathbf N}$ of the fractional Laplacian?

Recall that the eigenvalues and eigenfunctions solve \begin{align*} \left(-\frac{d^2}{d x^2}\right)^{\alpha / 2} \varphi(x)=\lambda \varphi(x), \quad x \in D, \end{align*} where $\varphi \in L^2(D)$ is extended to $\mathbf{R}$ by 0, and that the fractional Laplacian is given by \begin{align*} \left(-\frac{d^2}{d x^2}\right)^{\alpha / 2} f(x)=c_\alpha \mathrm{pv} \int_{-\infty}^{\infty} \frac{f(x)-f(y)}{|x-y|^{1+\alpha}} d y, \quad x \in \mathbf{R} \end{align*} where \begin{align*} c_\alpha=\frac{2^\alpha \Gamma\left(\frac{1+\alpha}{2}\right)}{\sqrt{\pi}\left|\Gamma\left(-\frac{\alpha}{2}\right)\right|} \end{align*}

It is known that the sequence satisfies $0<\lambda_1<\lambda_2 \leq \lambda_3 \leq \ldots$ and an estimate (but not their exact value) is given in https://arxiv.org/pdf/1012.1133.pdf.

Question 2: Can we also plot $\lambda_n +in$ in the complex plane?

---

Note: there are other questions on Mathematica SE that focus on plotting eigenfunctions of the Fractional Laplacian in multi-d, but I couldn't find anything about eigenvalues in 1-d.

Answer 1

We can compute eigenvalues with using Haar wavelets collocation method as follows (see also our code here)

c[n_, a_] := 2^(a) Gamma[(a + n)/2]/(Pi^(n/2) Abs[Gamma[-a/2]]);
lap[n_, a_, x_, u_] := 
  c[n, a] Integrate[(u[x] - u[y])/Abs[x - y]^(n + a), y];


x0 = -1; x1 = 1; nn = 40; dx = (x1 - x0)/(nn); xl = 
 Table[x0 + l*dx, {l, 0, nn}]; tcol = 
 Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}];
Psi[x_] = 
 Table[WaveletPhi[HaarWavelet[], (x - xl[[i]])/dx], {i, nn}]; var = 
 Array[v, {nn}]; u[t_] := var . Psi[t];
s = 1/2;

int = Table[
    Table[NIntegrate[(Psi[x][[j]] - Psi[y][[j]])/(x - y)^(1 + s), {y, 
       x0, x}, Method -> "PrincipalValue", Exclusions -> tcol], {j, 
      nn}], {x, tcol}] + 
   Table[Table[
     NIntegrate[(Psi[x][[j]] - Psi[y][[j]])/(y - x)^(1 + s), {y, x, 
       x1}, Method -> "PrincipalValue", Exclusions -> tcol], {j, 
      nn}], {x, tcol}];

int0 = Table[
   Table[Psi[x][[j]] NIntegrate[
      1/(x - y)^(1 + s), {y, -Infinity, x0}], {j, nn}], {x, tcol}];

int1 = Table[
   Table[Psi[x][[j]] NIntegrate[
      1/(y - x)^(1 + s), {y, x1, Infinity}], {j, nn}], {x, tcol}];

lp = c[1, s] (int + int0 + int1);
eq = Table[var . lp[[i]] == 0, {i, nn}]; {vec, matrix} = 
 CoefficientArrays[eq, var];

l = Eigenvalues[-matrix // N] // Reverse

Out[]= {-0.968028, -1.59898, -2.02729, -2.38745, -2.69731, \
-2.97784, -3.23278, -3.47001, -3.69072, -3.89884, -4.09479, -4.28086, \
-4.45717, -4.62515, -4.78483, -4.93714, -5.08206, -5.22023, -5.3516, \
-5.47662, -5.59523, -5.70774, -5.81411, -5.91456, -6.00903, -6.09769, \
-6.18049, -6.25753, -6.32878, -6.39432, -6.45411, -6.50821, -6.55659, \
-6.59928, -6.63627, -6.66757, -6.69317, -6.71309, -6.72732, -6.73585} 

These results we can compare with asymptotic formula from the paper

lambda = Table[(n Pi/2 - (2 - s) Pi/8)^(s) // N, {n, nn}]

Out[]= {0.990832, 1.59767, 2.0306, 2.38624, 2.69535, 2.9725, \
3.22591, 3.46083, 3.68078, 3.8883, 4.0853, 4.27323, 4.45324, 4.62624, \
4.793, 4.95416, 5.11023, 5.26168, 5.40889, 5.5522, 5.6919, 5.82825, \
5.96148, 6.0918, 6.21939, 6.34442, 6.46703, 6.58736, 6.70552, \
6.82165, 6.93582, 7.04815, 7.15872, 7.2676, 7.37488, 7.48062, \
7.58488, 7.68773, 7.78922, 7.88941}

In one plot they look like

ListPlot[{Abs[l], lambda}, PlotRange -> {0, lambda[[nn]]}, 
 PlotLegends -> {"Haar wavelets", "Asymptotic formula"}]

Therefore about 15 eigenvalues computed with error of about 1% on the grid with 40 collocation points. First six modes can be computed with Eigenvectors as

g = Eigenvectors[-matrix // N, -6];

Do[y[i] = 
   Interpolation[
    Join[{{-1, 0}}, 
     Table[{tcol[[j]], Re[g[[i]][[j]]]}, {j, nn}], {{1, 0}}]];, {i, 6}]

Table[Plot[y[i][x], {x, -1, 1}, Frame -> True], {i, 6}] // Reverse  

Question 2. Complex plot in a case of nn=80 and $\alpha =0.5$

ComplexListPlot[{Table[Abs[l[[n]]] + I n, {n, Length[l]}], 
  Table[lambda[[n]] + I n, {n, Length[l]}]}, AspectRatio -> 1, 
 Frame -> True, PlotLabel -> Row[{"\[Alpha] = ", .5}], 
 PlotLegends -> {"Haar wavelets", "Asymptotic formula"}]