5
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ What about parameters $a, b$ in equation definition? Is it typo that these parameters used in interval definition as end points, while in initial data there is unit interval? $\endgroup$ Apr 10, 2022 at 2:35

1 Answer 1

7
$\begingroup$

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

Figure 1

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 Figure 2

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. Figure 3

$\endgroup$
8
  • $\begingroup$ Thank you very much! This is great! Indeed, I had a typo in the question (the interval was supposed to be $(\alpha, \beta)$. What happens if it is e.g. (-1,1)? $\endgroup$
    – Riku
    Apr 12, 2022 at 13:48
  • $\begingroup$ @Riku See update to my answer with numerical solution on an arbitrary interval $x_0 < x < x_1$. $\endgroup$ Apr 12, 2022 at 17:25
  • $\begingroup$ Thank you very much! $\endgroup$
    – Riku
    Apr 12, 2022 at 20:54
  • $\begingroup$ @Riku You are welcome! See last update with 28 collocation points and maximal error about $1.18\times 10^-6$ $\endgroup$ Apr 13, 2022 at 2:28
  • $\begingroup$ Thanks for the update! That's great news. By the way, I tried to adapt your code to the hyperbolic case, but I couldn't make it work. I've posted the attempt to this post: mathematica.stackexchange.com/questions/266701/… If you have some time, could you please take a look? $\endgroup$
    – Riku
    Apr 13, 2022 at 20:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.