5
$\begingroup$

How can I use Mathematica to solve the following higher order and nonlinear fractional equation? $$\partial_t u(t,x) +\partial_x(u^k(t,x)\partial_x(-\Delta)^su(t,x)) = f(t,x) \quad t >0, \ x \in (\alpha,\beta),\\ u(t,x) = u_c \quad t \ge 0, \ x \in \mathbb{R} \setminus (\alpha,\beta), \\ (-\Delta)^su(t,x) = 0, \quad t \ge 0, \ x \in \{\alpha,\beta\}, \\ u(0,x) = u_0(x) \quad x \in (\alpha,\beta)$$ for a real number $k \ge 1$ and smooth functions u_0,u_c f. Here $(-\Delta)^s$ is the singular integral fractional Laplacian.

Note that the second-order and linear case was amazingly solved in Solve 1d fractional parabolic equations with Mathematica

Also note that a similar example (for only k = 1 and including an additional term) appears in this paper: https://arxiv.org/pdf/1611.00164.pdf

Additionally, I would also be interested in the case where $(-\Delta)^s$ is the spectral fractional Laplacian and the boundary conditions are replaced by $$\partial_x u(t,x) = u_c \quad t \ge 0, \ x \in \{\alpha,\beta\}, \\ u^k(t,x)\partial_x(-\Delta)^su(t,x) = 0, \quad t \ge 0, \ x \in \{\alpha,\beta\}, \\ $$

Improve this question
$\endgroup$

1 Answer 1

Reset to default
5
+50
$\begingroup$

We can solve this problem in two steps. First we compute fractional Laplacian $v=(-\Delta)^\frac{s}{2}u$, and then we solve fractional thin film equation $\partial_tu=\partial_x(u^k\partial_x v)$. Code for the step one is based on the colocation method with using Euler wavelets (see my previous answer here and here). Let define exact solution to test code

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[a_, x_] := (1 + x^2)^(-(1 - a)/2);
uc[x_, t_] := ue[s, x] (c0 kp (1/(c0 kp) - t))^(-1/kp); 
vc[x_, t_] := (c0 kp (1/(c0 kp) - t))^(-1/kp)
   2^s Gamma[(1 + s)/2]/Gamma[(1 - s)/2] (1 + x^2)^(-(1 + s)/2);
lape[a_, x_, t_] := (c0 kp (1/(c0 kp) - t))^(-1/kp)
    2^a Gamma[(1 + a)/2]/Gamma[(1 - a)/2] (1 + x^2)^(-(1 + a)/2);
f[x_, t_] := 
  1/Gamma[(1 - s)/2] (1 - c0 kp t)^(-1 - 1/kp) (1 + x^2)^(-(5/2) - s/
    2) (-c0 (1 + x^2)^(2 + s) Gamma[(1 - s)/2] + 
     2^s kp (-1 + s) (1 + s) (-1 + 
        c0 kp t) x^2 ((1 - c0 kp t)^(-1/kp) (1 + x^2)^(1/2 (-1 + s)))^
      kp Gamma[(1 + s)/2] - 
     2^s (1 + s) (3 + s) (-1 + 
        c0 kp t) x^2 ((1 - c0 kp t)^(-1/kp) (1 + x^2)^(1/2 (-1 + s)))^
      kp Gamma[(1 + s)/2] + 
     2^s (1 + s) (-1 + c0 kp t) (1 + 
        x^2) ((1 - c0 kp t)^(-1/kp) (1 + x^2)^(1/2 (-1 + s)))^
      kp Gamma[(1 + s)/2]);

Here kp=k and c0 some parameters defined below. On these lines we define wavelets and function to be solved

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 = 2; M0 = 3; x0 = 0; 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}]; xcol = 
 tcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 With[{k = k0, M = M0}, PsiE[k, M, t1]]; Intx1 = 
 With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Intx2 = Integrate[Intx1, t1]; 
intx1[y_] := Intx1 /. t1 -> (y - x0)/(x1 - x0); 
intx2[y_] := Intx2 /. t1 -> (y - x0)/(x1 - x0); 
Psi[y_] := Psijk /. t1 -> (y - x0)/(x1 - x0);

varu[t_] := Table[u[i][t], {i, nn}]; 
varv[t_] := Table[v[i][t], {i, nn}]; 
v[x_, t_] := varv[t] . intx2[x] + x v0[t] + v01[t]; 
v1[x_, t_] := varv[t] . intx1[x] + v0[t]; 
v2[x_, t_] := varv[t] . Psi[x]; 
u[x_, t_] := varu[t] . intx1[x] + u01[t];
u1[x_, t_] := varu[t] . Psi[x];

Code to compute fractional Laplacian

s = 9/10; kp = 3; c0 = 1/6; int = 
 Table[Table[
    NIntegrate[(intx1[tcol[[i]]][[j]] - 
        intx1[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}] // 
  Quiet;

int1 = Table[
    Table[NIntegrate[(intx1[tcol[[i]]][[j]] - 
         intx1[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}] // 
   Quiet;

int0 = Table[
    Table[NIntegrate[(intx1[tcol[[i]]][[
         j]])/(tcol[[i]] - y)^(1 + s), {y, -Infinity, x0}, 
      AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}] // 
   Quiet;

int2 = Table[
    Table[NIntegrate[(intx1[tcol[[i]]][[
         j]])/(y - tcol[[i]])^(1 + s), {y, x1, Infinity}, 
      AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}] // 
   Quiet;
lp = c[1, s] (int + int1 + int0 + int2) // Re;

System of equations, initial and boundary conditions

eq = Join[
   Table[D[varu[t], t] . lp[[i]] - D[v[xcol[[i]], t], t] == 0, {i, 
     nn}], Table[-D[u[xcol[[i]], t], t] + 
      u[xcol[[i]], t]^kp v2[xcol[[i]], t] + 
      kp u[xcol[[i]], t]^(kp - 1) u1[xcol[[i]], t] v1[xcol[[i]], t] - 
      f[xcol[[i]], t] == 0, {i, nn}]];
ic = Join[Table[u[tcol[[i]], 0] == uc[tcol[[i]], 0], {i, nn}], 
   Table[v[tcol[[i]], 0] == vc[tcol[[i]], 0], {i, nn}]];
bc = {v[x0, t] == vc[x0, t], v[x1, t] == vc[x1, t], 
  u[x0, t] == uc[x0, t]}; bct = {D[v[x0, t], t] == D[vc[x0, t], t], 
  D[v[x1, t], t] == D[vc[x1, t], t], 
  D[u[x0, t], t] == D[uc[x0, t], t]};
var = Join[{u01, v0, v01}, Table[v[i], {i, nn}], 
  Table[u[i], {i, nn}]]; ic0 = {u[x0, 0] == uc[x0, 0], 
  v[x0, 0] == vc[x0, 0], v[x1, 0] == vc[x1, 0]};

Finally we solve system of equations and visualize fractional Laplacian (exact and numeric solution) and error in collocation points

sol = NDSolve[{eq, ic, ic0, bct}, var, {t, 0, 1}, 
   Method -> {"EquationSimplification" -> "Residual"}];

{Plot3D[vc[x, t], {x, x0, x1}, {t, 0, 1}, ColorFunction -> "Rainbow", 
  Mesh -> None, PlotTheme -> "Marketing", AxesLabel -> Automatic, 
  PlotRange -> All], 
 Plot3D[v[x, t] /. sol[[1]], {x, x0, x1}, {t, 0, 1}, 
  ColorFunction -> "Rainbow", Mesh -> None, PlotTheme -> "Marketing", 
  AxesLabel -> Automatic, PlotRange -> All], 
 Plot[Evaluate[Table[vc[x, t] - v[x, t] /. sol, {x, tcol}]], {t, 0, 
   1}, PlotLegends -> tcol, AxesLabel -> Automatic, PlotRange -> All]}

Figure 1

It looks nice, but we need to compute $\partial_x v$, and it is why error increases drastically when we using $\partial_x v$ to compute u

lst = Table[{{x, t}, v1[x, t] /. sol[[1]]}, {x, 
    Join[{x0}, xcol, {x1}]}, {t, 0, 1, .01}];

vx = Interpolation[Flatten[lst, 1], InterpolationOrder -> 4];

U = NDSolveValue[{-D[w[x, t], t] + D[w[x, t]^kp vx[x, t], x] == 
    f[x, t], w[x, 0] == uc[x, 0], w[x0, t] == uc[x0, t], 
   w[x1, t] == uc[x1, t]}, w, {x, x0, x1}, {t, 0, 1}];

Visualization exact (left), numerical solution (center), and error on collocation points

{Plot3D[uc[x, t], {x, x0, x1}, {t, 0, 1}, ColorFunction -> "Rainbow", 
  Mesh -> None, PlotTheme -> "Marketing", AxesLabel -> Automatic, 
  PlotRange -> All], 
 Plot3D[U[x, t], {x, x0, x1}, {t, 0, 1}, ColorFunction -> "Rainbow", 
  Mesh -> None, PlotTheme -> "Marketing", AxesLabel -> Automatic, 
  PlotRange -> All], 
 Plot[Evaluate[Table[uc[x, t] - U[x, t], {x, tcol}]], {t, 0, 1}, 
  PlotLegends -> tcol, AxesLabel -> Automatic, PlotRange -> All]}

Figure 2 We also can ignore last step and compute u with using sol. In this case it looks even better then U

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

Figure 3

Update 1. Note, that in general case the thin film model describes waves. To reproduce wave solutions we can use as initial data ue and ve from the paper Finite difference methods for fractional Laplacians in the form

Clear["Global`*"]

ue[x_, k_, s_] := (1 - x^2)^(k + s/2);
lap[x_, k_, s_] := 
  Hypergeometric2F1[(1 + s)/2, -k, 1/2, x^2] 2^s Gamma[
    k + 1 + s/2] Gamma[(1 + s)/2]/(k! Gamma[1/2]);

The corresponding f can be evaluated as D[ue[x, k, s]^n D[lap[x, k, s], x], x], therefore, we have

f[x_, k_, s_, n_] := (1/(Sqrt[\[Pi]] k!))
   2^(2 + s)
     k n (k + s/2) (1 + s) x^2 (1 - x^2)^(-1 + k + s/2) ((1 - x^2)^(
     k + s/2))^(-1 + n)
     Gamma[1 + k + s/2] Gamma[(1 + s)/2] Hypergeometric2F1[1 - k, 
     1 + (1 + s)/2, 3/2, x^2] - (
   2^(1 + s) k (1 + s) ((1 - x^2)^(k + s/2))^
    n Gamma[1 + k + s/2] Gamma[(1 + s)/2] Hypergeometric2F1[1 - k, 
     1 + (1 + s)/2, 3/2, x^2])/(Sqrt[\[Pi]] k!) - (
   2^(3 + s) (1 - k) k (1 + s) (1 + (1 + s)/2) x^2 ((1 - x^2)^(
     k + s/2))^
    n Gamma[1 + k + s/2] Gamma[(1 + s)/2] Hypergeometric2F1[2 - k, 
     2 + (1 + s)/2, 5/2, x^2])/(3 Sqrt[\[Pi]] k!);

With this function we solve equation $\partial_tu=\partial_x(u^k\partial_x v)-f$ with using code

c[n_, a_] := a 2^(a - 1) Gamma[(a + n)/2]/(Pi^(n/2) Gamma[1 - a/2]);
res = {0 < a < 2};
lapc[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 = 2; M0 = 4; x0 = -1; 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}]; xcol = 
 tcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 With[{k = k0, M = M0}, PsiE[k, M, t1]]; Intx1 = 
 With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Intx2 = Integrate[Intx1, t1]; 
intx1[y_] := Intx1 /. t1 -> (y - x0)/(x1 - x0); 
intx2[y_] := Intx2 /. t1 -> (y - x0)/(x1 - x0); 
Psi[y_] := Psijk /. t1 -> (y - x0)/(x1 - x0);

varu[t_] := Table[u[i][t], {i, nn}]; 
varv[t_] := Table[v[i][t], {i, nn}]; 
v[x_, t_] := varv[t] . intx2[x] + x v0[t] + v01[t]; 
v1[x_, t_] := varv[t] . intx1[x] + v0[t]; 
v2[x_, t_] := varv[t] . Psi[x]; 
u[x_, t_] := varu[t] . intx1[x] + u01[t];
u1[x_, t_] := varu[t] . Psi[x];

s = 1/8; kp = 2; L = 2 (1 + s) c[1, s]; mL = 
 Gamma[2 + s/2] Gamma[5/2 + s/2]; uc[x_] := 1.5 ue[x, 1, s]; 
vc[x_] := lap[x, 1, s];

int = Table[
    Table[NIntegrate[(intx1[tcol[[i]]][[j]] - 
         intx1[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}] // 
   Quiet;

int1 = Table[
    Table[NIntegrate[(intx1[tcol[[i]]][[j]] - 
         intx1[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}] // 
   Quiet;

int0 = Table[
    Table[NIntegrate[(intx1[tcol[[i]]][[
         j]])/(tcol[[i]] - y)^(1 + s), {y, -Infinity, x0}, 
      AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}] // 
   Quiet;

int2 = Table[
    Table[NIntegrate[(intx1[tcol[[i]]][[
         j]])/(y - tcol[[i]])^(1 + s), {y, x1, Infinity}, 
      AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}] // 
   Quiet;
(*intb0=Table[NIntegrate[uc[y,0]/(tcol[[i]]-y)^(1+s),{y,-Infinity,x0}]\
,{i,nn}]//Quiet;
intb1=Table[NIntegrate[uc[y,0]/(y-tcol[[i]])^(1+s),{y,x1,Infinity}],{\
i,nn}];*)
lp = c[1, s] (int + int1 + int0 + int2) // Re;
(*eq1=v[x,t]==(-\[CapitalDelta])^(s/2)u[x,t];
eq2=-D[u[x,t],t]+u[x,t]^k D[v[x,t],x,x]+k u[x,t]^(k-1) \
D[u[x,t],x]D[v[x,t],x]+L x D[u[x,t],x]+L u[x,t];*)

eq = Join[
   Table[D[varu[t], t] . lp[[i]] - D[v[xcol[[i]], t], t] == 0, {i, 
     nn}], Table[-D[u[xcol[[i]], t], t] + 
      u[xcol[[i]], t]^kp v2[xcol[[i]], t] + 
      kp u[xcol[[i]], t]^(kp - 1) u1[xcol[[i]], t] v1[xcol[[i]], t] - 
      f[xcol[[i]], 1, s, kp] == 0, {i, nn}]];
ic = Join[Table[u[tcol[[i]], 0] == uc[tcol[[i]]], {i, nn}], 
   Table[v[tcol[[i]], 0] == vc[tcol[[i]]], {i, nn}]];
bc = {v[x0, t] == vc[x0, t], v[x1, t] == vc[x1, t], 
  u[x0, t] == uc[x0, t]}; bct = {D[v[x0, t], t] == 0, 
  D[v[x1, t], t] == 0, D[u[x0, t], t] == 0};
var = Join[{u01, v0, v01}, Table[v[i], {i, nn}], 
  Table[u[i], {i, nn}]]; ic0 = {u[x0, 0] == 0, v[x0, 0] == vc[x0], 
  v[x1, 0] == vc[x1]};

sol = NDSolve[{eq, ic, ic0, bct}, var, {t, 0, 1}, 
   Method -> {"EquationSimplification" -> "Residual"}];

Visualization numerical solution with very clear waves for u (left) and v (right)

{Plot3D[u[x, t] /. sol[[1]], {x, x0, x1}, {t, 0, 1}, 
  ColorFunction -> "Rainbow", Mesh -> None, PlotTheme -> "Marketing", 
  AxesLabel -> Automatic, PlotRange -> All, PlotPoints -> 50], 
 Plot3D[v[x, t] /. sol[[1]], {x, x0, x1}, {t, 0, 1}, 
  ColorFunction -> "Rainbow", Mesh -> None, PlotTheme -> "Marketing", 
  AxesLabel -> Automatic, PlotRange -> All, PlotPoints -> 50]}

Figure 4

Update 2. We can exclude step with fractional Laplacian computation and compute numerical solution in the first case 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[a_, x_] := (1 + x^2)^(-(1 - a)/2);
uc[x_, t_] := ue[s, x] (c0 kp (1/(c0 kp) - t))^(-1/kp); 
vc[x_, t_] := (c0 kp (1/(c0 kp) - t))^(-1/kp)
   2^s Gamma[(1 + s)/2]/Gamma[(1 - s)/2] (1 + x^2)^(-(1 + s)/2);
lape[a_, x_, t_] := (c0 kp (1/(c0 kp) - t))^(-1/kp)
    2^a Gamma[(1 + a)/2]/Gamma[(1 - a)/2] (1 + x^2)^(-(1 + a)/2);
f[x_, t_] := 
  1/Gamma[(1 - s)/2] (1 - c0 kp t)^(-1 - 1/kp) (1 + x^2)^(-(5/2) - s/
    2) (-c0 (1 + x^2)^(2 + s) Gamma[(1 - s)/2] + 
     2^s kp (-1 + s) (1 + s) (-1 + 
        c0 kp t) x^2 ((1 - c0 kp t)^(-1/kp) (1 + x^2)^(1/2 (-1 + s)))^
      kp Gamma[(1 + s)/2] - 
     2^s (1 + s) (3 + s) (-1 + 
        c0 kp t) x^2 ((1 - c0 kp t)^(-1/kp) (1 + x^2)^(1/2 (-1 + s)))^
      kp Gamma[(1 + s)/2] + 
     2^s (1 + s) (-1 + c0 kp t) (1 + 
        x^2) ((1 - c0 kp t)^(-1/kp) (1 + x^2)^(1/2 (-1 + s)))^
      kp Gamma[(1 + s)/2]);
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 = 2; M0 = 4; x0 = 0; 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}]; xcol = 
 tcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 With[{k = k0, M = M0}, PsiE[k, M, t1]]; Intx1 = 
 With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Intx2 = Integrate[Intx1, t1]; 
intx1[y_] := Intx1 /. t1 -> (y - x0)/(x1 - x0); 
intx2[y_] := Intx2 /. t1 -> (y - x0)/(x1 - x0); 
Psi[y_] := Psijk /. t1 -> (y - x0)/(x1 - x0);

varu[t_] := Table[u[i][t], {i, nn}];
u[x_, t_] := varu[t] . intx2[x] + u01[t] + u0[t] x;
u1[x_, t_] := varu[t] . intx1[x] + u0[t];
u2[x_, t_] := varu[t] . Psi[x];


s = 9/10; kp = 3; c0 = 1/6;
intu1 = Table[
    Table[NIntegrate[(intx1[tcol[[i]]][[j]] - 
         intx1[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}] // 
   Quiet;

int1u1 = Table[
    Table[NIntegrate[(intx1[tcol[[i]]][[j]] - 
         intx1[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}] // 
   Quiet;

int0u1 = Table[
    Table[NIntegrate[(intx1[tcol[[i]]][[
         j]])/(tcol[[i]] - y)^(1 + s), {y, -Infinity, x0}, 
      AccuracyGoal -> 10, PrecisionGoal -> 8], {j, nn}], {i, nn}] // 
   Quiet;

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

int1u2 = 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}] // 
   Quiet;

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

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

lpu1 = c[1, s] (intu1 + int1u1 + int0u1 + int2u1) // Re;(*v1*)
lpu2 = c[1, s] (intu2 + int1u2 + int0u2 + int2u2) // Re;(*v2*)


eq = Table[-D[u[xcol[[i]], t], t] + 
     u[xcol[[i]], t]^kp varu[t] . lpu2[[i]] + 
     kp u[xcol[[i]], t]^(kp - 1) u1[xcol[[i]], t] varu[t] . 
       lpu1[[i]] - f[xcol[[i]], t] == 0, {i, nn}];
ic = Table[u[tcol[[i]], 0] == uc[tcol[[i]], 0], {i, nn}];
bc = {u[x0, t] == uc[x0, t], 
  u[x1, t] == uc[x1, t]}; bct = {D[u[x1, t], t] == D[uc[x1, t], t], 
  D[u[x0, t], t] == D[uc[x0, t], t]};
var = Join[{u01, u0}, 
  Table[u[i], {i, nn}]]; ic0 = {u[x0, 0] == uc[x0, 0], 
  u[x1, 0] == uc[x1, 0]};

sol = NDSolve[{eq, ic, ic0, bct}, var, {t, 0, 1}, 
   Method -> {"EquationSimplification" -> "Residual"}];

Visualization exact and numerical solution, and error in the collocation points

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

Figure 5

Improve this answer
$\endgroup$
10
  • $\begingroup$ Thank you so much! How did you choose vc in this? That is, how is it related to the originally given uc? $\endgroup$
    – Riku
    Apr 21, 2022 at 18:07
  • $\begingroup$ @Riku Fractional Laplacian for ue[a_, x_] := (1 + x^2)^(-(1 - a)/2) has been computed in the paper Numerical Methods for the Fractional Laplacian: A Finite Difference-Quadrature Approach by Yanghong Huang and Adam Oberman. I used this ue to construct the time dependent separable solution uc[x, t] and corresponding vc[x,t] for fractional thin film equation. $\endgroup$
    – Alex Trounev
    Apr 22, 2022 at 1:26
  • $\begingroup$ Thank you. I'm trying to use (1 + x^2)^(-(1 - a)/2) as initial data of the equation and 2^s Gamma[(1 + s)/2]/Gamma[(1 - s)/2] (1 + x^2)^(-(1 + s)/2) as corresponding initial data for (-\Delta)^s u with external boundary data = 0 and f = 0, but I get as a result that the solution does not move; I'm probably making a mistake, but I don't see where. Could you maybe separate in the data of your question initial conditions from boundary conditions? $\endgroup$
    – Riku
    Apr 22, 2022 at 6:46
  • $\begingroup$ @Riku In the section System of equations, initial and boundary conditions there are eq (equations), ic (initial conditions), bc (boundary conditions), bct (differentiated boundary conditions), ic0 (initial conditions for bct). For usage fractional Laplacian of ue see my answer on mathematica.stackexchange.com/questions/266357/… $\endgroup$
    – Alex Trounev
    Apr 22, 2022 at 9:52
  • $\begingroup$ I see. It's very strange though: I tried to initialize the same initial condition ue but with zero boundary data and (mistakenly) got a stationary profile for all times 0 < t < 1. For completeness, could you please add a simulation where the initial data is just ue[a_, x_] := (1 + x^2)^(-(1 - a)/2) (and the corresponding ve[x_]:=2^s Gamma[(1 + s)/2]/Gamma[(1 - s)/2] (1 + x^2)^(-(1 + s)/2) and zero boundary condition? $\endgroup$
    – Riku
    Apr 22, 2022 at 22:17

Your Answer

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

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