0
$\begingroup$

I'm numerically integrating the following system,

f[r_] = 0.904965 - (1 + 12.7028/r^2) (1 - 2/r)
dϕdλ = 1/r[λ]^2;
eqc = r'[λ] == Sqrt[f[r[λ]]/(dϕdλ^2)];

with the numerical routine,

sc = NDSolve[{eqc, r[0] == 4.85744 + .003},r, {λ, 0., 40*Pi}, AccuracyGoal -> 10, PrecisionGoal -> 20, Method -> {"StiffnessSwitching"}]

However, I need to keep the integration when the solution reaches the radius r = 4.85744 or r = 11.3301, which are the values where the derivative changes the signal.

Technically, what I need is: When r[λ] = 4.85744 or r[λ] = 11.3301 keep the integration with r'[λ]==-Sqrt[f[r[λ]]/(dϕdλ^2)].

In my code, the integration just stops at this point. I would like to keep integrating with opposite velocity when the solution reaches these two radii, is it possible?

Thanks in advance.

$\endgroup$
  • $\begingroup$ There was a little mistake but I've fixed and edited the post, I hope I made it clear @ChrisK. $\endgroup$ – Herr Schrödinger May 14 at 20:14
  • 1
    $\begingroup$ Thanks for the edit! $\endgroup$ – Chris K May 14 at 20:15
  • $\begingroup$ There is a solution with reflection from point r=11.3301, but there is no solution with reflection from point r=4.85744. $\endgroup$ – Alex Trounev May 14 at 23:25
  • $\begingroup$ @AlexTrounev That’s exactly the solution that I want to integrate! $\endgroup$ – Herr Schrödinger May 14 at 23:43
2
$\begingroup$

It is necessary to square the equation eqc and differentiate it once by lambda. Then we can use WhenEvent[r[\[Lambda]] == 4.85744, r'[\[Lambda]] -> -r'[\[Lambda]]].

F[r_] := 2 (0.904965` - (1 + 12.7028`/r^2) (1 - 2/r)) r + ((
     25.4056` (1 - 2/r))/r^3 - (2 (1 + 12.7028`/r^2))/r^2) r^2
mu[r_] := Sqrt[0.904965 - (1 + 12.7028/r^2) (1 - 2/r)]*r;
eqc = 2*r''[\[Lambda]] == F[r[\[Lambda]]]; r0 = 4.86;
sc = NDSolve[{eqc, r[0] == r0, r'[0] == mu[r0], 
   WhenEvent[r[\[Lambda]] == 4.85744, 
    r'[\[Lambda]] -> -r'[\[Lambda]]]}, r, {\[Lambda], 0., 80*Pi}, 
  Method -> {"StiffnessSwitching"}]

Plot[Evaluate[r[\[Lambda]] /. sc], {\[Lambda], 0., 80*Pi}, 
 PlotRange -> {0, 12}, AxesLabel -> {"\[Lambda]", "r"}]

fig1

$\endgroup$

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.