# Solve hyperbolic fractional problem in 1d

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?

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