# Solve 1d fractional parabolic equations with Mathematica

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

• 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? Apr 10, 2022 at 2:35

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.

• 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)?
– Riku
Apr 12, 2022 at 13:48
• @Riku See update to my answer with numerical solution on an arbitrary interval $x_0 < x < x_1$. Apr 12, 2022 at 17:25
• Thank you very much!
– Riku
Apr 12, 2022 at 20:54
• @Riku You are welcome! See last update with 28 collocation points and maximal error about $1.18\times 10^-6$ Apr 13, 2022 at 2:28
• 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?
– Riku
Apr 13, 2022 at 20:05