Solve 1d fractional parabolic equations with Mathematica

Question

How can we use Mathematica to solve the following fractional parabolic problem? $$\partial_t u(t,x) +a(-\Delta)^su(t,x) + bu(t,x)= f(t,x) \quad t >0, \ x \in (\alpha,\beta),\\ u(t,x) = u_c(t,x)\quad t \ge 0, \ x \in \mathbb{R} \setminus (\alpha,\beta), \\ u(0,x) = u_0(x) \quad x \in (\alpha,\beta)$$ for $a,b\ge0$ and smooth functions u_c, u_0, f. Here $(-\Delta)^s$ is the singular integral fractional Laplacian.

The stationary case of this problem was brilliantly solved in a related post: Solve 1d fractional equation with Mathematica

Answer 1

The numerical algorithm is not so differ from that we described on this page. We can use collocation method and Euler wavelets as well as follows

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];

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;
s = 1/2; d = 1; b = 1; ue[a_, x_] := (1 + x^2)^(-(1 - a)/2); 
uc[x_, t_] := ue[s, x] Exp[-t];
lape[a_, x_, t_] := 
  Exp[-t] 2^a Gamma[(1 + a)/2]/Gamma[(1 - a)/2] (1 + x^2)^(-(1 + a)/2);


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, 0]/(tcol[[i]] - y)^(1 + s), {y, -Infinity, 0}], {i, nn}];

intb1 = Table[
      NIntegrate[
    uc[y, 0]/(y - tcol[[i]])^(1 + s), {y, 1, Infinity}], {i, 
        nn}];



var[t_] := Table[v[i][t], {i, nn}]; ic = 
 Table[u[tcol[[i]], 0] == uc[tcol[[i]], 0], {i, nn}];
u[x_, t_] := var[t] . Psi[x]; lp = 
 d c[1, s] (int + int1 + int0 + int2);

f[x_, t_] := lape[s, x, t] + b uc[x, t] - uc[x, t]; eq = 
 Table[D[u[tcol[[i]], t], t] + var[t] . lp[[i]] - 
    d c[1, s] Exp[-t] (intb0[[i]] + intb1[[i]]) + 
    b u[tcol[[i]], t] - (f[tcol[[i]], t]) == 0, {i, nn}];

sol = NDSolve[{eq, ic}, Table[v[i], {i, nn}], {t, 0, 1}];

Visualization numerical solution and difference exact and numerical solutions in collocation points. Note that maximal absolute error is about $1.17\times 10^{-6}$ for 16 collocation points.

{Plot3D[u[x, t] /. sol[[1]], {x, 0, 1}, {t, 0, 1}, 
  ColorFunction -> "Rainbow", Mesh -> None, PlotTheme -> "Marketing", 
  AxesLabel -> Automatic], 
 Plot[Evaluate[Table[uc[x, t] - u[x, t] /. sol, {x, tcol}]], {t, 0, 
   1}, PlotLegends -> tcol, AxesLabel -> Automatic]}

In a case of arbitrary interval $x_0 < x < x_1$ the code can be modified as follows (in this example $x_0=-2, x_1=1$)

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];

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; x0 = -2; x1 = 1; nn = 
 Total[With[{k = k0, M = M0}, 
   Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]]; dx = (x1 - 
    x0)/(nn); xl = Table[x0 + 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 - x0)/(x1 - x0);
s = 1/2; d = 1; b = 1; ue[a_, x_] := (1 + x^2)^(-(1 - a)/2);
uc[x_, t_] := ue[s, x] Exp[-t];
lape[a_, x_, t_] := 
  Exp[-t] 2^a Gamma[(1 + a)/2]/Gamma[(1 - a)/2] (1 + x^2)^(-(1 + a)/2);


int = Table[
   Table[NIntegrate[(Psi[tcol[[i]]][[j]] - 
        Psi[y][[j]])/(tcol[[i]] - y)^(1 + s), {y, x0, 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]], x1}, 
     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, x0}, AccuracyGoal -> 10, 
     PrecisionGoal -> 8], {j, nn}], {i, nn}];

int2 = Table[
   Table[NIntegrate[(Psi[tcol[[i]]][[j]])/(y - tcol[[i]])^(1 + s), {y,
       x1, Infinity}, AccuracyGoal -> 10, PrecisionGoal -> 8], {j, 
     nn}], {i, nn}];
intb0 = Table[
   NIntegrate[
    uc[y, 0]/(tcol[[i]] - y)^(1 + s), {y, -Infinity, x0}], {i, nn}];

intb1 = Table[
   NIntegrate[
    uc[y, 0]/(y - tcol[[i]])^(1 + s), {y, x1, Infinity}], {i, nn}];



var[t_] := Table[v[i][t], {i, nn}]; ic = 
 Table[u[tcol[[i]], 0] == uc[tcol[[i]], 0], {i, nn}];
u[x_, t_] := var[t] . Psi[x]; lp = 
 d c[1, s] (int + int1 + int0 + int2);

f[x_, t_] := lape[s, x, t] + b uc[x, t] - uc[x, t]; eq = 
 Table[D[u[tcol[[i]], t], t] + var[t] . lp[[i]] - 
    d c[1, s] Exp[-t] (intb0[[i]] + intb1[[i]]) + 
    b u[tcol[[i]], t] - (f[tcol[[i]], t]) == 0, {i, nn}];

sol = NDSolve[{eq, ic}, Table[v[i], {i, nn}], {t, 0, 1}];

Visualization of numerical solution and error in collocation points

Note, that in this case the maximal absolute error is about $3.14\times 10^{-4}$, and this is much higher then in the case of unit interval. To decrease error we can increase number of collocation points, for example from 16 to 28 with k0=3, M0=7. In this case the maximal absolute error is about $1.18\times 10^{-6}$, that practically is same as in previous example with unit interval.