3
$\begingroup$

I have a differential equation given by 

-4*D[S[u,v],u,v]==V[u,v]*S[u,v]

with boundary conditions 

S[0,v]==Exp[-(v-vc)^2/(2*sigma^2)]
S[u,0]==Exp[-(-vc)^2/(2*sigma^2)]

V[u,v]=(1-a/r-b*r^2)*((l*(1+l))/r^2+(a/r^2-2*b*r)/r)

where r should be determined from 

v-u=2*Integrate[1/(1-a/r-b*r^2),r]

(NOTE: l is a nonzero positive integer, and a and b are two nonzero positive numbers. "It is possible to use numbers instead of a, b, and l through calculations") Indeed, I know that I should use discretization to solve this equation because I solved a simpler case with b=0. As a point, I used 

First[NDSolve[
    {-4*D[S[u, v], u, v] == V[u, v]* S[u, v],S[u,0]==Exp[-(-vc)^2/2*sigma^2)],
     S[0,v]==Exp[-(v-vc)^2/(2*sigma^2)]},S, {u, 0, 400}, {v, 0, 400},
     Method -> {"MethodOfLines","SpatialDiscretization" -> 
     {"TensorProductGrid","MaxStepSize" -> 0.1}}, AccuracyGoal -> 10]]

to solve the simpler case. But there is a problem here. r cannot be calculated straightforwardly from v-u=2*Integrate[1/(1-a/r-b*r^2),r].

One can use a=1, l=1, b=1/16, vc=10, and sigma=3 for numerical analysis.

As the final remark, maybe it is possible to use Runge-Kutta method for the integration, interpolates the function at each step and finally finds r as a function of u and v. But I do not know how I can do it :(

I will be really thankful if someone help.

Improve this question
$\endgroup$
  • 1
    $\begingroup$ If this equation is to be solved numerically, please provide values for all constants. $\endgroup$ – bbgodfrey Jul 10 '18 at 2:27
  • $\begingroup$ @bbgodfrey a=1, l=1, b=1/16, vc=10, and sigma=3. $\endgroup$ – Mehrab Jul 10 '18 at 14:56
  • $\begingroup$ @bbgodfrey I checked that answer and found out it is wrong (see comment below it). If you have any ideas, please let me know. $\endgroup$ – Mehrab Jul 10 '18 at 20:09
  • $\begingroup$ v - u = 2*Integrate[1/(1-a/r-b*r^2), r] yields v - u = 2/(b r) for small v - u, which does not seem credible to me. Perhaps, you mean v - u = c + 2*Integrate[1/(1-a/r-b*r^2), r]. If so, what is c? Once this is resolved, solving the problem should be straightforward. $\endgroup$ – bbgodfrey Jul 10 '18 at 23:03
  • $\begingroup$ @bbgodfrey No, there is no c 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 for r[x] for all x because NDSolveValue cannot solve this equation. This is the solution: we should write a code to solve this equation with small steps from -x0 to x0, not this big interval {x, -x0, x0}. Did you get it? $\endgroup$ – Mehrab Jul 10 '18 at 23:39
5
$\begingroup$

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]

enter image description here

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

enter image description here

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]

enter image description here

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

enter image description here

Improve this answer
$\endgroup$
  • $\begingroup$ bbgodfrey, thank you for the solution of the problem. But how to get the data in Fig. 1 of this article using this solution? $\endgroup$ – Alex Trounev Jul 13 '18 at 4:42
  • $\begingroup$ @AlexTrounev Thank you for suggesting that I add the third plot to my answer. Note that the parameters used to obtain Figure 1 do not appear to match those contained in the question. $\endgroup$ – bbgodfrey Jul 14 '18 at 1:01
  • $\begingroup$ Thank you, your solution of the problem is very close to what we are discussing referring to the article. But there are a couple of mistakes in the article, one of which you pointed out, and I pointed to the other. The question arises, should we continue to refer to the article or stop at your solution? $\endgroup$ – Alex Trounev Jul 14 '18 at 3:00
  • $\begingroup$ @AlexTrounev I believe that we have answered the question but look forward to hearing from Mehrab. $\endgroup$ – bbgodfrey Jul 14 '18 at 3:50
  • $\begingroup$ @bbgodfrey Dear Godfrey, it seems to me that your answer to this question is correct mathematically. Thank you very much for spending time and solving mine problem. But let me help if you want to reproduce Fig. 1 of the paper. I was working on your code to reproduce it but I couldn't yet. However, I know how to reproduce it. The potential barrier must be real and positive in the tortoise coordinate between re and rc. Indeed, this problem should be solved for this range [re,rc], because it is our real accessible spacetime... Please see next comment... $\endgroup$ – Mehrab Jul 15 '18 at 0:50
3
$\begingroup$ 
    a = 1; b = 1/16; x0 = 800; r = 
 NDSolveValue[{D[r0[x], x] == (1 - a/r0[x] - b*r0[x]^2)/2, 
   r0[-x0] == 1 + I*10^-10}, r0, {x, -x0, x0}]
vc = 10; sigma = 3; l = 1;
S0[u_, v_] := Exp[-(v - vc)^2/(2*sigma^2)]
eq = {-4*D[S[u, v], u, v] == V[Re[r[v - u]]]*S[u, v]};
bc = {S[0, v] == S0[0, v], S[u, 0] == S0[u, 0]};
V[x_] := (1 - a/x - b*x^2)*((l*(1 + l))/x^2 + (a/x^2 - 2*b*x)/x)
sol = NDSolveValue[{eq, bc}, S, {u, 0, x0}, {v, 0, x0}]
{Plot3D[Re[sol[u, v]], {u, 0, x0}, {v, 0, x0}, PlotPoints -> 50, 
  Mesh -> None, PlotRange -> All, AxesLabel -> {"u", "v"}], 
 Plot[Re[r[x]], {x, -x0, -x0 + 10}, AxesLabel -> {"v-u", "r"}, 
  PlotStyle -> Orange, PlotRange -> All], 
 Plot[Re[r[x]], {x, -x0, x0}, AxesLabel -> {"v-u", "r"}, 
  PlotStyle -> Orange, PlotRange -> All]}

fig2

After numerous discussions with the author of the topic, I decided to add a solution to the problem on the basis of my numerical model, but with the data that Godfrey used. As can be seen from the data given, there is a coincidence and a difference. At first, the maximum of the function $V(x)$ is shifted to the left with respect to 0, and in the Godfrey solution the maximum is shifted to the right. The author of the theme asserts that the maximum should be shifted to the left. And he sent me the appropriate decision, with which I agree. Secondly, there are differences in solution, which is due to the first difference.

a = 1; b = 1/16; x0 = 100;
R1 = -4.42863948675507037483099540169716097781583927443899886866276463\
30182560928665325495919425006839287038680064401788392`100.;
R2 = 1.078377745621778233051751805397040016199774219084253403438456893\
3353720068855339101075303006206231499526837887379443`100.;
R3 = 3.350261741133292141779243596300120961616065055354745465224307739\
6828840859809986394844122000633055539153226514408949`100.;
ndvar = -16`100 ((Log[r - R1] R1)/(-16 + 3 R1^2) + (
     Log[r - R2] R2)/(-16 + 3 R2^2) + (Log[R3 - r] R3)/(-16 + 3 R3^2));
r00 = r /. 
    FindRoot[ndvar == 0, {r, 11/10`100}, AccuracyGoal -> Infinity, 
     MaxIterations -> 10000, PrecisionGoal -> 100, 
     WorkingPrecision -> 100]; // Quiet
r = NDSolveValue[{D[r0[x], x] == (1 - a/r0[x] - b*r0[x]^2)/2, 
   r0[0] == r00}, r0, {x, -x0, x0}]
vc = 10; sigma = 3; l = 1;
S0[u_, v_] := Exp[-(v - vc)^2/(2*sigma^2)]
eq = {-4*D[S[u, v], u, v] == V[Re[r[v - u]]]*S[u, v]};
bc = {S[0, v] == S0[0, v], S[u, 0] == S0[u, 0]};
V[x_] := (1 - a/x - b*x^2)*((l*(1 + l))/x^2 + (a/x^2 - 2*b*x)/x)
sol = NDSolveValue[{eq, bc}, S, {u, 0, x0}, {v, 0, x0}]
{Plot3D[Re[sol[u, v]], {u, 0, x0}, {v, 0, x0}, PlotPoints -> 50, 
  Mesh -> None, PlotRange -> All, AxesLabel -> {"u", "v"}], 
 Plot[V[Re[r[x]]], {x, -50, 50}, AxesLabel -> {"v-u", "V"}, 
  PlotStyle -> Orange, PlotRange -> All], 
 Plot[Re[r[x]], {x, -20, 60}, AxesLabel -> {"v-u", "r"}, 
  PlotStyle -> Orange, PlotRange -> All], 
 LogLogPlot[Abs[sol[50, v]], {v, 5, x0}]}

fig2

Improve this answer
$\endgroup$
  • $\begingroup$ Good work. But did you check the values of r[x] for different x's? They are the same! Also, you have applied the wrong boundary conditions. Thank you for your time. $\endgroup$ – Mehrab Jul 10 '18 at 20:05
  • $\begingroup$ This answer 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 for r[x] for all x because NDSolveValue cannot solve this equation. This is the solution: we should write a code to solve this equation with small steps from -x0 to x0, not this big interval {x, -x0, x0}. But I do not know how I can do this. $\endgroup$ – Mehrab Jul 10 '18 at 23:46
  • $\begingroup$ I showed the working code, and you can play with the parameters, substitute any initial values and look for solutions. Good luck! $\endgroup$ – Alex Trounev Jul 11 '18 at 2:39
  • $\begingroup$ Mehrab, I did not see your answer with the data for the parameters when I was preparing my code. I made some improvement to cover this case and published a new code and picture.. $\endgroup$ – Alex Trounev Jul 11 '18 at 5:23
  • $\begingroup$ Dear friend, I really tried to explain and open this problem and tell it in an easy way. But obviously I could not, and I know that it was my weakness that could not explain the problem well. I am really thankful for the generous help and spending time. As the final point, I want you to have a look at pages 3 and 4 of arxiv.org/pdf/gr-qc/9902010v2.pdf. Maybe it explains this problem easier. Best wishes. $\endgroup$ – Mehrab Jul 11 '18 at 7:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.