0
$\begingroup$

I am trying to solve these two equations for various values of B(magnetic field): $$\text{xb}'(l)=\frac{\sqrt{\text{g11}\left(r_0\right) \left(-\text{gtt}\left(r_0\right)\right)}}{\text{g11}(\text{rb}(l)) \sqrt{-\text{gtt}(\text{rb}(l))}};$$ $$\text{rb}'(l)=\frac{\sqrt{\text{g11}\left(r_0\right) \text{gtt}\left(r_0\right)-\text{g11}(\text{rb}(l)) \text{gtt}(\text{rb}(l))}}{\sqrt{\text{g11}(\text{rb}(l)) \text{grr}(\text{rb}(l)) (-\text{gtt}(\text{rb}(l)))}};$$ where:

f[r_] = 1 - rh^4/r^4 - (2 B^2)/(3 r^4) Log[r/rh];

q[r_] = 1 - (2 B^2)/(3 r^4) Log[r];

h[r_] = 1 + B^2/(3 r^4) Log[r];

gtt[r_] = -r^2 f[r];

g11[r_] = r^2 h[r];

grr[r_] = 1/(r^2 f[r]);

For B=1, when I used NDSolve for the second equation, I only get half of the answer:

xb1[l_] = Sqrt[-gtt[r0] g11[r0]]/(Sqrt[-gtt[rb[l]]] g11[rb[l]]);

rb1[l_] = 
  Sqrt[-gtt[rb[l]] g11[rb[l]] + gtt[r0] g11[r0]]/
   Sqrt[-gtt[rb[l]] g11[rb[l]] grr[rb[l]]];
B = 1;

rh = 1;

r0 = 1.1;
rs = NDSolve[{y'[l] == rb1[l] /. rb -> y, y[0] == 1.11}, y, {l, -2, 2}]

NDSolve::mxst: Maximum number of 99124 steps reached at the point l == -0.101218.

{{y -> InterpolatingFunction[{{-0.101218, 2.}}, <>]}}
Output for rb[l] Is there any algorithm to get the full solution? (If possible, an analytical one?)

Improve this question
$\endgroup$
6
  • 2
    $\begingroup$ For y<1 rb1is complex! $\endgroup$
    – Ulrich Neumann
    Feb 22, 2022 at 15:08
  • $\begingroup$ @UlrichNeumann: I think the transition is actually at y = 1.1, which is the value of r0 that's been programmed in. The square root in the denominator of rb1 certainly vanishes when rb[l] == r0, and I expect that it switches from positive to negative there. $\endgroup$
    – Michael Seifert
    Feb 22, 2022 at 18:39
  • $\begingroup$ @UlrichNeumann I am trying to reproduce results from this paper: sciencedirect.com/science/article/pii/S0370269322000831 On page 3 of the paper, they have written "rb (l) and xb (l) are the static solutions obtained integrating the above two equations", the results of which are further used in calculations. Is there any error in interpreting their statement? $\endgroup$
    – codebpr
    Feb 23, 2022 at 5:03
  • $\begingroup$ @MichaelSeifert yes that's true but when I use the condition rb[0]=r0, I get a straight line instead of the above plot. $\endgroup$
    – codebpr
    Feb 23, 2022 at 5:06
  • $\begingroup$ @codebpr Sorry, I don't know the details in this paper and didn't find a hint how the Plots y<1.1are created. $\endgroup$
    – Ulrich Neumann
    Feb 23, 2022 at 8:08

2 Answers 2

Reset to default
2
$\begingroup$

In the linked paper, the equation of motion for $\text{rb}$ has a plus-or-minus sign in front of it; equivalently, it can be written as a second-degree, first-order ODE: $$ (\text{rb}'(l))^2=\frac{\text{g11}\left(r_0\right) \text{gtt}\left(r_0\right)-\text{g11}(\text{rb}(l)) \text{gtt}(\text{rb}(l))}{\text{g11}(\text{rb}(l)) \text{grr}(\text{rb}(l)) (-\text{gtt}(\text{rb}(l)))}. $$ If we want to find for a solution with a continuous value of $\text{rb}'$ and $\text{rb}''$, as is shown in Figure 3 of the linked paper, we can differentiate both sides of the equation to yield a higher-order ODE: $$ 2 (\text{rb}'(l)) (\text{rb}''(l)) = {\text{a huge mess, but one that} \choose \text{ Mathematica can calculate} } $$ To ensure that this is equivalent to our original EOM, we must enforce the original lower-order EOM at $l = 0$ as an initial condition: $$ \text{rb}'(0)= \pm \sqrt{\frac{\text{g11}\left(r_0\right) \text{gtt}\left(r_0\right)-\text{g11}(\text{rb}(0)) \text{gtt}(\text{rb}(0))}{\text{g11}(\text{rb}(0)) \text{grr}(\text{rb}(0)) (-\text{gtt}(\text{rb}(0)))}}. $$

To implement this, run the setup code and then the following:

neweom = D[(y'[l])^2 == (rb1[l])^2 /. rb -> y, l];
newrs = NDSolve[{neweom, y'[0] == rb1[0] /. rb -> y, y[0] == 1.11}, 
  y, {l, -2, 2}]
Plot[y[l] /. newrs, {l, -2, 2}]

enter image description here

Improve this answer
$\endgroup$
4
  • $\begingroup$ It may also be possible to get Mathematica to use different numerical solving algorithms to yield this solution automatically, but I was unable to find a good way to do this. In particular, I had hoped that Method -> {"EquationSimplification" -> "Residual", "IndexReduction" -> {Automatic, "IndexGoal" -> 0}} seemed would do the trick, but it yielded a solution with a discontinuous value of $\text{rb}'$ at $l = 0$. $\endgroup$
    – Michael Seifert
    Feb 23, 2022 at 15:48
  • $\begingroup$ Thank you for this illuminating Answer. If I put y[0]=r0, I get a complex infinity. So I used y[0]=1.10001 and it works perfectly. I have to further use it for solving other equations. Is there a way to use this solution for rb[l] from NDSolve as an input to a different equation? $\endgroup$
    – codebpr
    Feb 24, 2022 at 5:26
  • $\begingroup$ @codebpr: I believe so, yes, though I'd have to play around with the code to figure out how to do it. That said, if you want to use it as a "background" for an ODE over the same region, the simplest way to do it is simply to solve for this function and the new one as a simultaneous system of ODEs. $\endgroup$
    – Michael Seifert
    Feb 24, 2022 at 13:37
  • $\begingroup$ Many thanks for the suggestion. I will try to apply it in the code. $\endgroup$
    – codebpr
    Feb 24, 2022 at 13:57
0
$\begingroup$

To long for a comment:

"Graphical solution" of the ode

StreamPlot[ {1, rb1[l] /. rb -> y /. y[l] -> y} , {l, -2, 2}, {y, 0.5,2}, FrameLabel -> {l, y}, GridLines -> {None, {1.1}}]

enter image description here

shows solutions only for y>1.1 (see comment @MichaelSeifert )

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for the answer. Actually, in the assumptions, it is given that r0=1.1 which is the tip of the plot -> ars.els-cdn.com/content/image/… $\endgroup$
    – codebpr
    Feb 23, 2022 at 11:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.