# 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? Commented Dec 22, 2021 at 16:14
• @Liuvv Where did you get this model? Commented Dec 22, 2021 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 Commented Dec 23, 2021 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. Commented Dec 23, 2021 at 5:36
• You should add these new info into the body of the question by clicking the Edit button. Commented Dec 23, 2021 at 10:32

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 :) . Commented Dec 25, 2021 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! Commented Dec 25, 2021 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. Commented Dec 25, 2021 at 13:02
• @xzczd Thank you, it works. Commented Dec 25, 2021 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 Commented Dec 30, 2021 at 11:17

# 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? Commented Dec 25, 2021 at 16:40
• @alex I just fixed OP's code. Commented Dec 26, 2021 at 4:40
• Thank you very much. Then we can try to improve this code. Commented Dec 26, 2021 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] Commented Dec 26, 2021 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. Commented Dec 26, 2021 at 14:05