Before beginning, it is important to observe that Eq. (2.10) in the article cited by the OP is incorrect. Because r0 < 0
, as stated earlier in the article, the argument of Log[r/r0 - 1]
is strictly negative, which is nonphysical. In fact, the term should be Log[1 - r/r0]
. This observation is important, because one of the stated goals of the question is to determine an analogous result for the parameters in the question.
Begin by performing the integration in the question relating v - u
to r
.
int = 2 Integrate[1/(1 - a/r - b*r^2) /. {a -> 1, b -> 1/16}, r] // Normal
(* -32 ((Log[r - Root[16 - 16 #1 + #1^3 &, 1]] Root[16 - 16 #1 + #1^3 &, 1])
/(-16 + 3 Root[16 - 16 #1 + #1^3 &, 1]^2)
+ (Log[r - Root[16 - 16 #1 + #1^3 &, 2]] Root[16 - 16 #1 + #1^3 &, 2])
/(-16 + 3 Root[16 - 16 #1 + #1^3 &, 2]^2)
+ (Log[r - Root[16 - 16 #1 + #1^3 &, 3]] Root[16 - 16 #1 + #1^3 &, 3])
/(-16 + 3 Root[16 - 16 #1 + #1^3 &, 3]^2)) *)
plus a constant of integration to be determined. Singularities occur at
Cases[int, _Root, 3];
%//N
(* {-4.42864, 1.07838, 3.35026} *)
which readily are identified as the quantities {r0, re, rc}
appearing in the article. According to the article, r0 == -(re + rc)
, which also is true here.
FullSimplify[Total@Rest@%% == -First@%%]
(* True *)
Next, transform int
to the form of the corrected Eq. (2.10).
dvar = int /. {Log[r - z_] /; z > 0 -> Log[r/z - 1],
Log[r - z_] /; z < 0 -> Log[-r/z + 1]} /.
Log[-1 + r/Root[16 - 16 #1 + #1^3 &, 3]] -> Log[1 - r/Root[16 - 16 #1 + #1^3 &, 3]]
(* -32 ((Log[1 - r/Root[16 - 16 #1 + #1^3 &, 1]] Root[16 - 16 #1 + #1^3 &, 1])
/(-16 + 3 Root[16 - 16 #1 + #1^3 &, 1]^2)
+ (Log[-1 + r/Root[16 - 16 #1 + #1^3 &, 2]] Root[16 - 16 #1 + #1^3 &, 2])
/(-16 + 3 Root[16 - 16 #1 + #1^3 &, 2]^2)
+ (Log[1 - r/Root[16 - 16 #1 + #1^3 &, 3]] Root[16 - 16 #1 + #1^3 &, 3])
/(-16 + 3 Root[16 - 16 #1 + #1^3 &, 3]^2)) *)
plus constants terms proportional to Root[16 - 16 #1 + #1^3 &, -]
and to I Pi
that can be absorbed into the constant of integration, which we now set equal to zero to obtain the expression matching the right side of the corrected Eq. (2.10). A plot of dvar
follows.
ParametricPlot[{dvar, r}, {r, -10, 10}, AspectRatio -> 1/GoldenRatio,
AxesLabel -> {"u - v", r}, ImageSize -> Large,
LabelStyle -> {Black, Bold, Medium}, PlotRange -> All]

Next, compute the potential as a function of v - u
.
V = Simplify[(1 - a/r - b*r^2)*((l*(1 + l))/r^2 + (a/r^2 - 2*b*r)/r)
/. {a -> 1, b -> 1/16} /. l -> 1];
ndvar = N@dvar;
fV[z_?NumericQ] := If[Abs[z] < 50, Re[V /. FindRoot[ndvar == z, {r, 11/10}]], 0];
Plot[fV[x], {x, -50, 50}, PlotRange -> All, AxesLabel -> {"v - u", "V"},
ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

The exponential decrease in V
for large Abs[v - u]
justifies setting it to zero for Abs[v - u] > 50
. Doing so save much computing time. Finally, compute and plot the solution of the PDE.
sol = NDSolveValue[{-4*D[S[u, v], u, v] == fV[v - u]*S[u, v],
S[u, 0] == Exp[-(-vc)^2/(2*sigma^2)],
S[0, v] == Exp[-(v - vc)^2/(2*sigma^2)]} /. {vc -> 10, sigma -> 3},
S, {u, 0, 100}, {v, 0, 100}, MaxStepSize -> .1];
Plot3D[sol[u, v], {u, 0, 100}, {v, 0, 100}, AxesLabel -> {u, v, S},
ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}, PlotRange -> All,
PlotPoints -> 200]

A slice through v - u
space at constant u
, plotted logarithmically shows oscillations that otherwise are hard to see.
LogLogPlot[Abs[sol[50, v]], {v, 5, 100}, ImageSize -> Large,
LabelStyle -> {Black, Bold, Medium}]

a=1
,l=1
,b=1/16
,vc=10
, andsigma=3
. $\endgroup$ – Mehrab Jul 10 '18 at 14:56v - u = 2*Integrate[1/(1-a/r-b*r^2), r]
yieldsv - u = 2/(b r)
for smallv - u
, which does not seem credible to me. Perhaps, you meanv - u = c + 2*Integrate[1/(1-a/r-b*r^2), r]
. If so, what isc
? Once this is resolved, solving the problem should be straightforward. $\endgroup$ – bbgodfrey Jul 10 '18 at 23:03c
beside the integral. I am sure about this. I know the method that should be used to solve the problem, but I do not know how I can apply it. The answer is given below is good but it has a problem.r=NDSolveValue[{D[r0[x], x] == (1 - a/r0[x] - b*r0[x]^2)/2, r0[-x0] == 1}, r0, {x, -x0, x0}]
gives an unique value forr[x]
for allx
becauseNDSolveValue
cannot solve this equation. This is the solution: we should write a code to solve this equation with small steps from-x0
tox0
, not this big interval{x, -x0, x0}
. Did you get it? $\endgroup$ – Mehrab Jul 10 '18 at 23:39