Solve hyperbolic fractional problem in 1d

Question

Inspired by the solution of a fractional parabolic problem presented in https://mathematica.stackexchange.com/posts/266475/, I wonder if we can also use Mathematica to solve the following fractional hyperbolic problem $$\partial_t^2 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)\\ \partial_tu(0,x) = u_1(x), \quad x \in (\alpha,\beta)$$ for $a,b\ge0$ and smooth functions u_c, u_0, u_1, f. Here $(-\Delta)^s$ is the singular integral fractional Laplacian.

I think that one should use

eq = 
 Table[D[u[tcol[[i]], t], {t,2}] + 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}];

to incorporate the second derivative and

ic = {Table[u[tcol[[i]], 0] == uc[tcol[[i]], 0], {i, nn}], Table[Derivative[x0,x1]u[tcol[[i]], 0] == u1[tcol[[i]], 0], {i, nn}]};

to incorporate the second initial data, but this produces many errors that I cannot decifer.

What is the correct way to proceed?

Answer 1

We can use last code from my answer here. We modified exact solution and add initial condition ic1, then we have

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 = 7; 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] Cos[t];
lape[a_, x_, t_] := 
  Cos[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}] // 
   Quiet;

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

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

int2 = Table[
    Table[NIntegrate[(Psi[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}];



var[t_] := Table[v[i][t], {i, nn}]; ic = 
 Table[u[tcol[[i]], 0] == uc[tcol[[i]], 0], {i, nn}]; ic1 = 
 Table[var'[0] . Psi[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, t] + var[t] . lp[[i]] - 
    d c[1, s] Cos[t] (intb0[[i]] + intb1[[i]]) + 
    b u[tcol[[i]], t] - (f[tcol[[i]], t]) == 0, {i, nn}];


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

Visualization of numerical solution and error in collocation points

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