# Plot eigenvalues of fractional Laplacian in 1D

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.

• Can you put all your equations into plain text MMA code for a specific sample problem? Sep 22, 2022 at 11:25
• @MarcoB Indeed, half of the core problem of this question is exactly how to implement the 1d fractional Laplacian into Mathematica code
– Riku
Sep 22, 2022 at 11:43
• You may want to make that more explicit in the question then. It currently seems to focus completely on the plotting, whereas the problem seems to start well before that. Sep 22, 2022 at 12:32
• There are several questions about fractional Laplacian on this Forum. For example, to compute $\lambda_i$ we can use code from my answer on mathematica.stackexchange.com/questions/266357/… Sep 22, 2022 at 15:27
• @AlexTrounev Thanks! How can that code be adapted to compute the $\lambda_i$?
– Riku
Sep 22, 2022 at 19:01

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"}]


• Amazing! Thank you so much! Just a question: how can I modify the line of code ListPlot[{Abs[l], lambda}, PlotRange -> {0, lambda[[nn]]}, PlotLegends -> {"Haar wavelets", "Asymptotic formula"}] to plot $\lambda_n + in$ in the complex plane instead?
– Riku
Sep 23, 2022 at 15:51
• @Riku See update to mu answer. Sep 23, 2022 at 16:25
• Thank you! I'm a bit confused though. What is plotted on the two axes of the complex plot? And what do the different colors of the points represent?
– Riku
Sep 23, 2022 at 22:35
• @Riku Actually, we plot $\lambda _n+ i n$, therefore in horizontal line we have $\lambda_n$ (numeric and asymptotic), and in vertical axis $i n$. The plot has been updated for a case nn=80. Sep 24, 2022 at 3:38