# Implement fractional Laplacian

What is a way to implement the Fractional Laplacian with Mathematica?

How can we apply such implementation to numerically solve the problem

$$(-\Delta)^su = 1 \text{ in } B_1(0), \\ u = 0 \text{ in } \mathbb{R}^2 \setminus B_1(0),$$ where $B_1(0)$ is the unit disk centered at the origin?

• Use momentum space. That is use the Fourier transform Restate boundary conditions as momentum space constraints. At end switch back to coordinate space, I.e. use inverse Fourier transform – user14199 Sep 26 '17 at 19:24
• Jay, what is your application for the Fractional Laplacian? – user21 Sep 26 '17 at 23:15

I can provide an approach by finite elements and an application of the functional calculus of selfadjoint operators.

# Background

The spectrum of the Laplacian on a bounded domain $$\varOmega$$ with sufficiently smooth boundary subject to homogeneous Dirichlet boundary conditions is discrete. By the spectral theorem, there are eigenfunctions $$e_i \in H^1_0(\varOmega)$$ and $$\lambda_i \in \mathbb{R}$$ with $$-\Delta \, e_i = \lambda_i \, e_i$$ and $$0 < \lambda_1 \leq \lambda_2 \leq \dotsc$$ and $$\int_\varOmega e_i(x) \, e_j(x) \, \operatorname{d} x = \delta_{ij}.$$ By functional calculus, the solution $$u \in H^{2s}(\varOmega) \cap H^{s}_0(\varOmega)$$ to $$(-\Delta)^s \,u = f$$ for $$f \in L^2(\varOmega)$$ can be represented by $$u = \sum_{i=1}^\infty \left( \lambda_i^{-s} \, e_i \cdot \int_\varOmega f(x) \, e_i(x) \, \operatorname{d} x \right).$$ Weyl's law states that the eigenvalues behave as $$\lambda_i \sim i^\frac{2}{\operatorname{dim}(\varOmega)}$$ for $$i \to \infty$$. (You can get an idea of the law by performing a $$\sin$$-Fourier transform on the eigenvalue equation $$- \Delta v = \lambda \, v$$ for functions on the square $$[0, 2\pi] \times [0, 2\pi]$$.)

Thus, it is meaningful to use a truncated expansion in order to approximate $$u$$.

Our aim is now to use finite elements to approximate the eigenfunctions $$e_i$$ and the eigenvalues $$\lambda_i$$.

# Implementation

First, we create a nice domain. Disks are boring so I use the following starfish:

R = DiscretizeRegion[
BoundaryMeshRegion[
Map[
t \[Function] (2 + Cos[5 t])/3 {Cos[t], Sin[t]},
Most@Subdivide[0., 2. Pi, 2000]],
Line[Partition[Range[2000], 2, 1, 1]]
],
MaxCellMeasure -> 0.001,
MeshQualityGoal -> "Maximal"
]


Next, we have to set up the finite element method. This is somewhat a mess but we cannot apply NDSolve directly; we need access to mass matrix and stiffness matrix of the system. Note that I use Sin[2 x] + Cos[x + 3 y] as right hand side $$f$$.

 Needs["NDSolveFEM"]
(*Initialization of Finite Element Method*)

Rdiscr = ToElementMesh[R, "MeshOrder" -> 1];
vd = NDSolveVariableData[{"DependentVariables", "Space"} -> {{u}, {x, y}}];
sd = NDSolveSolutionData[{"Space"} -> {Rdiscr}];
cdata = InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {{-IdentityMatrix[2]}},
"MassCoefficients" -> {{1}},
"LoadCoefficients" -> {{Sin[2 x] + Cos[x + 3 y]}}
];
bcdata = InitializeBoundaryConditions[vd, sd, {{DirichletCondition[u[x, y] == 0., True]}}];
mdata = InitializePDEMethodData[vd, sd];

(*Discretization*)
dpde = DiscretizePDE[cdata, mdata, sd];
dbc = DiscretizeBoundaryConditions[bcdata, mdata, sd];
{load, stiffness, damping, mass} = dpde["All"];


Having mass and stiffness matrix, we have to reduce them to the interior degrees of freedom. For FEM of order 1, the Dirichlet boundary conditions are deployed into the stiffness matrix by replacing rows that correspond to boundary vertices by the respective row of an identity matrix. Moreover, we exploit that Mathematica sorts boundary vertices in front.

(*Finding interior degrees of freedom*)

intplist = Min[UpperTriangularize[stiffness, 1]["NonzeroPositions"][[All, 1]]] ;;;

stiffness2 = stiffness[[intplist, intplist]];
mass2 = mass[[intplist, intplist]];


We need eigenvalues and eigenvectors. Since calculating them is expensive, we restrict our attention to the most relevant 5 percent.

(*Spectral decomposition and functional calculus*)
s = 3/4;
reducedmodeldimension = Floor[Length[stiffness2] 0.05];
{Λ, U} = Reverse /@ Eigensystem[{stiffness2, mass2},
-reducedmodeldimension,
Method -> {"Arnoldi", "MaxIterations" -> 4000
}]; // AbsoluteTiming
U = Map[u \[Function] u/Sqrt[u.mass2.u], U];


The last step is necessary to ensure that the row vectors of U (representing the eigenfunctions of $$-\Delta$$) form an $$L^2$$-orthonormal system (the $$L^2$$-inner product being represented by the reduced mass matrix mass2).

For example, we can draw the first 6 eigenfunction like this:

GraphicsGrid[
Partition[
Table[
eigenvec = ConstantArray[0, Dimensions[stiffness][[2]]];
eigenvec[[intplist]] = Flatten[U[[i]]];
eigenfun = ElementMeshInterpolation[{Rdiscr}, eigenvec];
Plot3D[eigenfun[x, y], {x, y} ∈ R],
{i, 1, 6}], 3],
ImageSize -> Large]


For the hidden operator L as defined in the comment below, we can easily optain its $$L^2$$-Moore-Penrose pseudoinverse Lpinv by

(*L=mass2.Transpose[Λ^s U].(U.mass2);*)
Lpinv = b \[Function] Transpose[U].(Λ^-s (U.b));


Next, we solve the reduced equations and write the result into the interior degrees of freedom of a zero vector. Finally, ElementMeshInterpolation provides us with an InterpolatingFunction that can be easily plotted.

solution = ConstantArray[0, Dimensions[stiffness][[2]]];
solfun = ElementMeshInterpolation[{Rdiscr}, solution];
Plot3D[solfun[x, y], {x, y} ∈ R]


As plausibility check: This is the result for $$(-\Delta)^s u = 1$$ (upon revisiting this post, I was puzzled by the solution being assymetric while OP asked for the solution of $$(-\Delta)^s u = 1$$ which ought to have the same symmetries as the domain):

# Discussion

This method is however very slow: It has complexity $$O(N^3)$$, where $$N$$ is the number of interior vertices. Moreover, this method is not very accurate.

Kernel based methods (e.g., convolution with fundamental solution) might be more accurate and might have complexity $$O(N^2)$$, but they might be limited to special domains where the fundamental solutions are known.

Speaking about special domains: For the disk $$B(0;1) \subset \mathbb{R}^2$$, the surface of the sphere $$S^2 \subset \mathbb{R}^3$$, the unit ball $$B(0;1) \subset \mathbb{R}^3$$, the standard tori $$\mathbb{T}^n = (S^1)^n$$, and all cube-like domains like $$Q = \coprod_{i=1}^n [a_i,b_i]$$, the eigenvectors and eigenvalues are well studied. E.g., using FFT will speed up the process tremendously for $$\mathbb{T}^n$$ and $$Q$$.

One could also circumvent the spectral decomposition step by directly assembling the Gagliardo bilinear form $$(u,v) \mapsto \pi^{-(2 s + n/2)}\frac{\Gamma(n/2+s)}{\Gamma(-s)}\int_\varOmega\int_\varOmega \frac{(u(x)-u(y))\,(v(x)-v(y))}{| x-y |^{n + 2 s}}\, \operatorname{d}x \, \operatorname{d}y$$ instead of the Laplacian stiffness matrix ($$n=2$$) (I am not exactly sure about the multiplicative constant in this formula; I got it from Hitchhiker's Guide to Fractional Sobolev Spaces, p. 8). This would require (partially singular) double integrals over $$O(N^2)$$ pairs of triangle elements. Moreover, that matrix is dense and needs $$O(N^3)$$ for factorization. Maybe one can apply FFT on a bounding box of the domain in order to construct a reasonable preconditioner for the conjugate gradient method.

• Thank you for your answer: it is very interesting. However, I'm not familiar with finite elements method (and actually I'm new to numerical analysis in general); would you mind adding some details on mass matrix, stiffness matrix, and the key ideas of your approach? – user52420 Sep 25 '17 at 13:31
• Thanks a lot. Could you point out a reference for the following statement? By functional calculus, the solution $u \in H^{s/2}_0(\varOmega)$ to $(-\Delta)^s \,u = f$ for $f \in L^2(\varOmega)$ can be represented by $$u = \sum_{i=1}^\infty \left( \lambda_i^{-s} \, e_i \cdot \int_\varOmega f(x) \, e_i(x) \, \operatorname{d} x \right).$$ Moreover, it is known that the eigenvalues behave as $\lambda_i \sim i^2$ for $i \to \infty$. – user52420 Sep 26 '17 at 10:24
• I have to confess that the assymptotics are not as good as state before. That is corrected now and I also learnt that this is called Weyl's law. – Henrik Schumacher Sep 26 '17 at 11:09
• Probably you wanted intplist without the Min – user21 Sep 27 '17 at 1:56
• @user21 Thank you for the comment. You see the three semicoli at the end of the line? The first two of them turn the minimum into a Span so that intplist acts indeed as a List when used in combination with Part. – Henrik Schumacher Sep 27 '17 at 5:50