Numerical solutions of active 1D wave equations

I would like to solve the following PDE with finite difference method. The PDE is from the following paper, (https://arxiv.org/pdf/1911.11823.pdf). I would like to implement the algorithm for the left panel of Figure 3 in this paper.

Below is a brief overview of the computational part of this paper.

$$-\frac{\partial ^{2}\Psi}{\partial t^{2}}+\frac{\partial ^2\Psi}{\partial {r^2_*}}-V_s(r_*)\Psi=0.$$.......................(1)

For this PDE, let's do the following coordinate transformation first,

$$u=t-r_*, v=t+r_*$$.

So we can get,

$$r_*=\frac{v-u}{2}$$,

$$r_*$$ is the tortoise coordinate. It can be calculated as following form,

$$dr_*=\frac{1}{f(r)}dr=\frac{1}{1-\frac{2M}{\sqrt{r^2+a^2}}}dr$$

And another function is $$V_s(r_*)$$, it is the effective potential, and we only know the following form,

$$V_s(r)=(1-\frac{2M}{\sqrt{r^2+a^2}})(\frac{l(l+1)}{r^2+a^2})$$

So $$V_s(r_*)$$ can be solved by using inverse function.

With the setup above, the equation (1) can be written as

$$4\frac{\partial ^{2}\psi(u,v)}{\partial u\partial v}=-V_s(r_*)\psi(u,v)$$.................... (2)

For this PDE, the appropriate discretization scheme is

$$\frac{\partial ^{2}\psi(u,v)}{\partial u\partial v}\to \frac{-\psi(u,h+v)-\psi(h+u,v)+\psi(h+u,h+v)+\psi(u,v)}{h^2}$$...............(3)

$$h$$ is the step of the grid. Then the equation (2) can be written as

$$4\frac{(-p(u,h+v)-p(h+u,v)+p(h+u,h+v)+p(u,v))}{h^2}+V_s\left(\frac{v-u}{2}\right) p(u,v)=0$$...........(4)

For this equation (4), Its initial conditions we know are

$$\psi(0,v)=e^{-\frac{(v-10)^2}{18}}, \psi(u,0)=0$$.

And I want to get the time evolution of $$\psi$$ and $$t$$.

Like the following picture, the values we used is $$M=0.5,a=1.01,l=1$$.

And the following is my try using Mathematica, but I failed. I would appreciate it if you could sort it out. And I'll fill you in on any details if you need. Thank you!

Clear["*"]
m = 1/2;
a = 1.01;
rs[r_] = Integrate[1/(1 - (2*m)/Sqrt[r^2 + a^2]), r];
l = 1;
V[r_] = (1 - (2*m)/Sqrt[r^2 + a^2]) (l (l + 1)/(r^2 + a^2));
Vs[r_] = V@InverseFunction[rs][r]; (*$$Vs[r_*]$$*)

{p[u + h, v + h] == p[u + h, v + h] + O[h]^3,
p[u + h, v] == p[u + h, v] + O[h]^3,
p[u, v + h] == p[u, v + h] + O[h]^3} // Normal

(* discretization *)
disc = Simplify[
Solve[Normal[{p[u + h, v + h] == p[u + h, v + h] + O[h]^3,
p[u + h, v] == p[u + h, v] + O[h]^3,
p[u, v + h] == p[u, v + h] + O[h]^3}], {D[p[u, v], u, v]}, {D[
p[u, v], u], D[p[u, v], v]}
]]

equForm = 4*D[p[u, v], u, v] + Vs[(v - u)/2] p[u, v] == 0 /. disc[[1]]

h = 1/2;
(*PDE Linear algebra*)
ecu[{u_, v_}] = equForm;

(* initial condition *)
p[0, v_] = Exp[-(v - 10)^2/(2*3*3)];
p[u_, 0] = 0;

(* Inside the solution range,the grid is divided *)(* the value of 17 \
in the following maybe is not right*)
coords = Flatten[
Table[{u, v}, {u, h, 17 - h, h/2}, {v, h, 17 - h, h/2}], 1];

(*PDE turns to Linear algebraic equations made up of many，many \
equations*)
ecus = ecu /@ coords;

(*An unknown quantity in a linear algebraic equation*)
vars = Union[Cases[ecus, p[_, _], \[Infinity]]];

(*Numerical solution of linear algebraic equations*)
sol = NSolve[ecus, vars];

(*I want to get the time evlution of \[Psi] and t, but I don't how to \
do it *)
LogPlot[Abs[p[0, v]], {v, 0, 17}, PlotRange -> All]


At the same time, I also share another way to calculate the $$\psi$$ and $$t$$ directly. The following is my try.

Clear["*"]
m = 1/2;
a = 1.01;
rs[r_] = Integrate[1/(1 - (2*m)/Sqrt[r^2 + a^2]), r];
l = 1;
V[r_] = (1 - (2 m)/Sqrt[r^2 + a^2]) (l (l + 1)/(r^2 + a^2));
Vs[r_] = V@InverseFunction[rs][r];
Plot[Vs[r], {r, -80, 80}, PlotRange -> All]

sol = NDSolveValue[{D[u[t, x], t, t] + Vs[x]*u[t, x] ==
D[u[t, x], x, x], u[0, x] == E^(-((x - 10)^2)/18),
Derivative[1, 0][u][0, x] == 0, u[t, -250] == 0, u[t, 250] == 0},
u, {t, 0, 500}, {x, -250, 250},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 50, "MinPoints" -> 50, "DifferenceOrder" -> 4}}](* Method\[Rule]{"ExplicitRungeKutta","StiffnessTest"\
\[Rule]True},MaxSteps\[Rule]1*^5,InterpolationOrder\[Rule]All *)

LogPlot[Abs[sol[t, 10]], {t, 0, 500}, PlotRange -> All]



However, campared with the figure of the paper, the accuracy of this method is not high. I think it's because the internal approach is not chosen properly. I tried using Runge Kutta method and the optimized discrete scheme provided by xzczd. This only requires changing the Method in NDSolve. It was for precision reasons that I chose to use the finite difference method. However, when using a coordinate transformation, the initial conditions will change.

• Mathematica has a framework for solving differential equations centered around NDSolve. Have you tried that out, instead of rolling your own solver? What led you to avoid the built in? Dec 22 '21 at 16:14
• @Liuvv Where did you get this model? Dec 22 '21 at 16:45
• Well, your @ won't work in this case. (As to the usage of @ you may want to read this: meta.stackexchange.com/q/43019/284701) Then, I'm in a hurry at the moment so can't go through the code, but 2 issues I can spot: 1. "I think I didn't put in boundary conditions… " No, after the coordinate transform the 2 i.c.s are enough to determine a solution even for FDM. 2. "I want to get the time evlution of \[Psi] and t… ", you need to transform from the $(u, v)$ coordinates back to the $(t, r_*)$ coordinates. BTW it this your classmate?: mathematica.stackexchange.com/q/260844/1871 Dec 23 '21 at 1:42
• “I tried to solve it using NDSolve with Runge Kutta method directly. ” You should not. The default ODE solver of NDSolve is quite robust, and the options for adjusting ODE solver should always be the last thing to touch. If the result isn't desired, try adjusting the options for spatial discretization first, here is an example: mathematica.stackexchange.com/a/196898/1871 "This is a method mentioned in a paper. " It's better to add the link of the paper to the question if it's convenient for you. Dec 23 '21 at 5:36
• You should add these new info into the body of the question by clicking the Edit button. Dec 23 '21 at 10:32

NDSolve-based Solution

tend = 500;
lb = -150; rb = -lb;
m = 1/2;
a = 1.01;
rs[r_] = Integrate[1/(1 - (2*m)/Sqrt[r^2 + a^2]), r];
l = 1;
V[r_] = (1 - (2 m)/Sqrt[r^2 + a^2]) (l (l + 1)/(r^2 + a^2));
Vs[r_] = V@InverseFunction[rs][r];
interVs = FunctionInterpolation[Vs[r], {r, lb, rb}]

mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}

abc = D[u[t, x], x] + direction D[u[t, x], t] ==
0 /. {{x -> lb, direction -> -1}, {x -> rb, direction -> 1}};

With[{Vs = interVs},
sol = NDSolveValue[{D[u[t, x], t, t] + Vs[x]*u[t, x] == D[u[t, x], x, x],
u[0, x] == E^(-((x - 10)^2)/18), Derivative[1, 0][u][0, x] == 0, abc},
u, {t, 0, tend}, {x, lb, rb}, Method -> mol[200, 2]]; // AbsoluteTiming]
(* {5.27039, Null} *)

LogPlot[sol[t, 10] // Abs, {t, 0, tend}, PlotRange -> All]


DensityPlot[sol[t, x], {t, 0, tend}, {x, lb, rb}, PlotRange -> All, PlotPoints -> 200,


Remark

1. Since the wave equation is solved in infinte space, I've set 1st order absorbing boundary condition (ABC) to relief reflection of the wave at the boundary. As shown in 2nd picture, though reflections still exist, it's not that large.

2. Further check shows that, even with naive Dirichlet boundary condition, sol[t, 10] // Abs remains almost the same at least for $$t\in[0,500]$$, $$x\in[-150,150]$$.

3. I've used FunctionInterpolation to rebuild Vs to avoid the inefficient InverseFunction.

4. The obtained solution isn't exactly the same as the one in the question, but the specific value of $$r_*$$ for FIG. 3 doesn't seem to be given in the paper. Also, the initial conditions for FIG. 3 aren't explicitly given in the paper either, so I'd like to stop here. Anyway, adjustions of spatial grid size, etc. show the solution above is probably reliable.

Quick Fix for the FDM Code of OP

As shown in the answers of Alex and mine, NDSolve is capable of handling the problem, but anyway, let me fix OP's FDM code. At least 4 issues here:

1. The most fatal one: Definition of coords is incorrect. The grids on $$u=0$$, $$v=0$$ are not covered.

2. The order of difference scheme you chose is $$1+1$$, which is rather low. (Notice the scheme in the paper is a $$2+2$$ order scheme. ) Since I'm aiming at a quick fix, I'll leave it alone and simply set a dense enough grid, but do notice if good performance is required, a better scheme is needed.

3. NSolve isn't bad, but for large linear system, LinearSolve is better.

4. Defintion of vars doesn't need to be that complicated.

Clear[h]
disc = Simplify[
Solve[Normal[{p[u + h, v + h] == p[u + h, v + h] + O[h]^3,
p[u + h, v] == p[u + h, v] + O[h]^3, p[u, v + h] == p[u, v + h] + O[h]^3}], {D[
p[u, v], u, v]}, {D[p[u, v], u], D[p[u, v], v]}]]

uend = 500/Sqrt[2];
vend = 500/Sqrt[2];
interVs2 = FunctionInterpolation[Vs[r], {r, -uend, uend}];

equForm = 4*D[p[u, v], u, v] + interVs2[(v - u)/2] p[u, v] == 0 /. disc[[1]]

h = 1/4;
ecu[{u_, v_}] = equForm;

p[0, v_] = N@Exp[-(v - 10)^2/(2*3*3)];
p[u_, 0] = 0.;

coords = Flatten[Table[{u, v}, {u, 0, uend - h, h}, {v, 0, vend - h, h}],
1]; // AbsoluteTiming

ecus = ecu /@ coords; // Quiet // AbsoluteTiming

vars = p @@@ ({h, h} + # & /@ coords); // AbsoluteTiming

{barray, marray} = CoefficientArrays[ecus, vars]; // AbsoluteTiming

sollst = LinearSolve[marray, -barray]; // AbsoluteTiming

solmat = Partition[sollst, uend/h // Round];

ArrayPlot[solmat, PlotRange -> All, DataReversed -> True,


solfunc = ListInterpolation[solmat, {{h, uend}, {h, vend}},
InterpolationOrder -> 1]; // AbsoluteTiming

With[{r = 37},
LogPlot[solfunc[t - r, t + r] // Abs, {t, r + h, 500/Sqrt[2] - r}, PlotRange -> All]]


Notice h = 1/4 may still be too coarse, but limited by the RAM of my laptop, I can't test further.

• It is a nice code (+1). Have you FDM code as well? Dec 25 '21 at 16:40
• @alex I just fixed OP's code. Dec 26 '21 at 4:40
• Thank you very much. Then we can try to improve this code. Dec 26 '21 at 4:59
• @xzczd It is very difficult to imagine that this problem will be solved in this way. Thank you for your help. I tried your code. After using interpolation, the speed of the program does become faster. However, in testing, I found that there is an error between the interpolated and unpolluted values. This error may differ in the presentation of the result. So, what do you think of this interpolation? Plot[{Va[x], Vs[x]}, {x, -100, 100}, PlotRange -> All] Dec 26 '21 at 13:36
• @Stone “I found that there is an error between the interpolated and unpolluted values. This error may differ in the presentation of the result” Add MaxRecursion -> 8 to definition of Va will make the discrepency negligible. Even with the original Va, I don't think the solution will be significantly influenced. Dec 26 '21 at 14:05

This problem can be solved with method of lines and with FDM as well. Using NDSolve we have

m = 1/2; int =
Integrate[1/(1 - (2*m)/Sqrt[r^2 + a^2]), r, Assumptions -> a > 0];
rs[x_] := int /. r -> x;
a = 1.01;
l = 1;
r[x_] := InverseFunction[rs][x];
V[r_] := (1 - (2*m)/Sqrt[r^2 + a^2]) (l (l + 1)/(r^2 + a^2));
Vs[x_] := V@InverseFunction[rs][x];


We can compare different definitions of potential

{Plot[r[x], {x, 0, 40}], Plot[rs[x], {x, 0, 10}],
Plot[V[r[x]], {x, 0, 30}], Plot[Vs[x], {x, 0, 30}]}


There are several options for method of lines to solve this problem. First variant is similar to FDM

eqs = {4*D[p[u, v], u, v] + Vs[(v - u)/2] p[u, v] == 0,
p[0, v] == Exp[-(v - 10)^2/(2*3*3)] Tanh[3 v], p[u, 0] == 0};
sol = NDSolveValue[eqs, p, {u, 0, 400}, {v, 0, 400},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 200, "MaxPoints" -> 1200,
"DifferenceOrder" -> 2}}]


Visualization in large and small scale

{LogPlot[Abs[sol[t - 37, t + 37]], {t, 37, 400 - 37}],
DensityPlot[Abs[sol[u, v]], {u, 0, 400}, {v, 0, 400},
PlotPoints -> 100, ColorFunction -> Hue, FrameLabel -> Automatic,
PlotRange -> All],
Plot3D[Abs[sol[u, v]], {u, 0, 100}, {v, 0, 100}, PlotPoints -> 100,
ColorFunction -> Hue, AxesLabel -> Automatic, PlotRange -> All,
Mesh -> None, Boxed -> False, PlotTheme -> "Scientific"]}


Second variant has high difference order, and less points

sol = NDSolveValue[eqs, p, {u, 0, 400}, {v, 0, 400},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 400, "MaxPoints" -> 800,
"DifferenceOrder" -> 8}}]


Third variant has same points as first one, but more high difference order

sol = NDSolveValue[eqs, p, {u, 0, 400}, {v, 0, 400},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 400, "MaxPoints" -> 1200,
"DifferenceOrder" -> 4}}]


Practically we have same pictures in 3 different options. To speedup computation we can use interpolation for inverse function as suggested by xzczd in the form

Va = FunctionInterpolation[Evaluate[Vs[x]], {x, -400, 400}] // Quiet

eqs = {4*D[p[u, v], u, v] + Va[(v - u)/2] p[u, v] == 0,
p[0, v] == Exp[-(v - 10)^2/(2*3*3)] Tanh[3 v], p[u, 0] == 0};
sol = NDSolveValue[eqs, p, {u, 0, 400}, {v, 0, 400},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 400, "MaxPoints" -> 1200,
"DifferenceOrder" -> 4}}] // AbsoluteTiming


It takes 4 second only. Now we can test FDM. We use algorithm described in the paper cited as follows

Clear["*"]
m = 1/2; int =
Integrate[1/(1 - (2*m)/Sqrt[r^2 + a^2]), r, Assumptions -> a > 0];
rs[x_] := int /. r -> x;
a = 1.01;
l = 1;
r[x_] := InverseFunction[rs][x];
V[r_] := (1 - (2*m)/Sqrt[r^2 + a^2]) (l (l + 1)/(r^2 + a^2));
Vs[x_] := V@InverseFunction[rs][x];
L = 100; Va =
FunctionInterpolation[Evaluate[Vs[x]], {x, -L, L}] // Quiet;
n = 1000; P = Array[p, {n, n}];
h = L/(n - 1); h2 = 1/8 h^2; eq =
Table[-P[[i + 1, j + 1]] + P[[i + 1, j]] + P[[i, j + 1]] -
P[[i, j]] -
h2 (Va[h (j - i - 1)/2] P[[i + 1, j]] +
Va[h (j - i + 1)/2] P[[i, j + 1]]) == 0, {i, 1, n - 1}, {j, 1,
n - 1}]; bc =
Join[Table[P[[1, j]] == Exp[-(h (j - 1) - 10)^2/(2*3*3)], {j, n}],
Table[P[[i, 1]] == 0, {i, 2, n}]];
eqs = Join[Flatten[eq], bc];
var = Flatten[Table[P[[i, j]], {i, n}, {j, n}]];
{b, m} = CoefficientArrays[eqs, var];

sol = LinearSolve[m, -b];
sol1 = Transpose[Partition[sol, n]];


We can compare FDM solution with NDSolve

eqs = {4*D[psi[u, v], u, v] + Va[(v - u)/2] psi[u, v] == 0,
psi[0, v] == Exp[-(v - 10)^2/(2*3*3)] Tanh[3 v], psi[u, 0] == 0};
sol2 = NDSolveValue[eqs, psi, {u, 0, 100}, {v, 0, 100},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 1000, "MaxPoints" -> 1000,
"DifferenceOrder" -> 4}}]


Visualization

{ArrayPlot[Abs[sol1], PlotRange -> All, ColorFunction -> Hue,
PlotLegends -> Automatic, DataReversed -> True, Frame -> False,
PlotLabel -> "FDM"],
DensityPlot[Abs[sol2[u, v]], {u, 0, 100}, {v, 0, 100},
PlotPoints -> 200, ColorFunction -> Hue, Frame -> False,
PlotRange -> All, PlotLegends -> Automatic, MaxRecursion -> 2,
PlotLabel -> "NDSolve"]}


Note that solutions look very similar, but nevertheless they have a small difference of $$4\times 10^{-3}$$

{ListPlot[{Table[sol1[[i, i]], {i, n}],
Table[sol2[x, x], {x, 0, L, h}]}, Frame -> True,
PlotLegends -> {"FDM", "NDSolve"}],
ListPlot[{Table[sol1[[i, i]], {i, n}] -
Table[sol2[x, x], {x, 0, L, h}]}, Frame -> True]}


Finally we put L=400; n=2000 in FDM code and compare FDM solution with NDSolve with options

sol2 = NDSolveValue[eqs, psi, {u, 0, 400}, {v, 0, 400},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 2000, "MaxPoints" -> 2000,
"DifferenceOrder" -> 4}}];


• You can use FunctionInterpolation to rebuild the Vs to speed up the calculation :) . Dec 25 '21 at 7:34
• @Alex Trounev Thank you for your thoughtful answer. From your answer, I learn a lot. And the question will be asked ia that why the initial condition include the tanh[3 v]? Secondly, it is quickly to calculate the Vs, but the sol is too slow. Can we have a way to make the programe more quick than now it is? At last, except the the method of Inverse function, can we have another to calculate theVs? Thank you! Dec 25 '21 at 12:46
• Tanh[3 v] added for consistency of boundary conditions. Actually it is not affected the solution and included to satisfy the NDSolve logic only. Dec 25 '21 at 13:02
• @xzczd Thank you, it works. Dec 25 '21 at 16:11
• @UlrichNeumann “Examples of "MethodOfLines" usually need two boundary conditions concerning the spatial coordinate v`” No… usually the number of b.c. is equal to the order of derivative in corresponding direction. Notice I'm not solving the same PDE as Alex. Related: math.stackexchange.com/q/450367/58219 Dec 30 '21 at 11:17