How can we use Mathematica to solve the following fractional elliptic problem? $$a(-\Delta)^su + bu= f(x) \text{ in } (0,1),\\ u = c \text{ in } \mathbb{R} \setminus (0,1), $$ for $a,b,c\ge0$ and a smooth function f. Here $(-\Delta)^s$ is the singular integral fractional Laplacian.
There is a related question Implement fractional Laplacian but it refers to a less general problem (a=1, b=0, $f\equiv 1$, c = 0) and uses a spectral definition of the fractional Laplacian.
This problem can be solved with using colocation method and Euler wavelets. We define fractional Laplacian as in the paper in the form
c[n_, a_] := a 2^(a - 1) Gamma[(a + n)/2]/(Pi^(n/2) Gamma[1 - a/2]);
lap[n_, a_, x_, u_] :=
c[n, a] Integrate[(u[x] - u[y])/Abs[x - y]^(n + a), y];
with $0 < a < 2$. Let consider equation
$$d (-\Delta)^\frac{a}{2}u(x)+b u(x)=f(x)$$ for $0<x<1$, and boundary conditions $$u=u_c(x), x\le0||x\ge 1$$ Our method is based on wavelets definition
UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 4; nn =
Total[With[{k = k0, M = M0},
Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]]; dx =
1/(nn); xl = Table[ l*dx, {l, 0, nn}]; tcol =
Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Psi[y_] := Psijk /. t1 -> y;
It means that in this particular case we use nn=16 colocation point. To compute numerical matrix for the fractional Laplacian we need first define parameters s=a, d, b and functions uc, f, for example,
s = 1/2;d=1;b=1; ue[a_, x_] := (1 + x^2)^(-(1 - a)/2);uc[x_] := ue[s, x];
lape[a_, x_] :=
2^a Gamma[(1 + a)/2]/Gamma[(1 - a)/2] (1 + x^2)^(-(1 + a)/2);
f[x_] := d lape[s, x] + b ue[s, x];
Then we compute
int =
Table[Table[
NIntegrate[(Psi[tcol[[i]]][[j]] -
Psi[y][[j]])/(tcol[[i]] - y)^(1 + s), {y, 0, tcol[[i]]},
Method -> "PrincipalValue", Exclusions -> tcol[[i]] - y == 0,
AccuracyGoal -> 6, PrecisionGoal -> 6], {j, nn}], {i, nn}];
int1 = Table[
Table[NIntegrate[(Psi[tcol[[i]]][[j]] -
Psi[y][[j]])/(y - tcol[[i]])^(1 + s), {y, tcol[[i]], 1},
Method -> "PrincipalValue", Exclusions -> y - tcol[[i]] == 0,
AccuracyGoal -> 6, PrecisionGoal -> 6], {j, nn}], {i, nn}];
int0 =
Table[Table[
NIntegrate[(Psi[tcol[[i]]][[
j]])/(tcol[[i]] - y)^(1 + s), {y, -Infinity, 0},
AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}];
int2 = Table[
Table[NIntegrate[(Psi[tcol[[i]]][[j]])/(y - tcol[[i]])^(1 + s), {y,
1, Infinity}, AccuracyGoal -> 10, PrecisionGoal -> 8], {j,
nn}], {i, nn}];
intb0 =
Table[NIntegrate[
uc[y]/(tcol[[i]] - y)^(1 + s), {y, -Infinity, 0}], {i, nn}];
intb1 = Table[
NIntegrate[uc[y]/(y - tcol[[i]])^(1 + s), {y, 1, Infinity}], {i,
nn}];
With this definitions we can prepare system of algebraic equations to be solved
var = Array[v, {nn}];
u[t_] := var . Psi[t]; lp = d c[1, s] (int + int1 + int0 + int2);
eq = Table[
var . lp[[i]] - d c[1, s] (intb0[[i]] + intb1[[i]]) +
b u[tcol[[i]]] - f[tcol[[i]]] == 0, {i, nn}];
This system can be solved with
sol = FindRoot[eq, Table[{var[[i]], 1/10}, {i, nn}]]
Visualization of numerical solution (red points), exact solution taken from the paper cited (blue line), and absolute error
{Show[Plot[{ue[s, t]}, {t, 0, 1}, PlotRange -> All,
PlotStyle -> Blue],
ListPlot[Table[{tcol[[i]], Re[u[tcol[[i]]] /. sol]}, {i, nn}],
PlotStyle -> Red]],
ListPlot[Table[Abs[ue[s, t] - Re[u[t] /. sol]], {t, tcol}],
PlotRange -> All, Filling -> Axis]}
Note, that maximal error is about $1.88\times 10^{-6}$ for 16 collocation points.
