3
$\begingroup$

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?

$\endgroup$

1 Answer 1

3
$\begingroup$

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

Figure 1

$\endgroup$

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.