# Solving a differential equation involving the square of the derivative

I want to solve the differential equation $$\left( \frac{dr}{d\lambda} \right)^2 = 1-\frac{L^2}{r^2} \left( 1-\frac{1}{r} \right)$$ where $L$ is some parameter. The behavior I am looking for is when $r$ starts "at infinity" at $\lambda=0$, reaches some minimum $r$, and increases out to infinity again. (I'm trying to describe the path that light takes in the presence of a Schwarzschild black hole.)

The issue I'm running into is when I use ParametricNDSolve, I can't tell Mathematica to smoothly transition from negative $dr/d\lambda$ to positive $dr/d\lambda$ after reaching the minimum $r$.

f[r_]:=1-1/r;
soln = ParametricNDSolve[{r'[\[Lambda]]^2 == 1 - L^2/r[\[Lambda]]^2 *f[r[\[Lambda]]], r[0] == 1000}, r, {\[Lambda], 0, 1000}, {L}]
Plot[r[30][\[Lambda]] /. soln, {\[Lambda], 0, 1000}]


In the above code, it plots the graph fine. However if I try to increase the range of $\lambda$ (say to 1200), it breaks down; the situation where I solve for $r'(\lambda)$ first and take the negative square root is identical. I'm not sure how to capture this extra information of transitioning from the negative to positive square root in the differential equation.

Sorry if this is a basic question, I'm rather new to Mathematica (and this site).

• I don't have time at the moment to write a complete answer, but I would look into using an EventLocator to "catch" the events where you flip from one branch to the other. Have a look here: reference.wolfram.com/language/tutorial/… Apr 13 '17 at 22:06
• The equation can be solved symbolically using λ == Integrate[1/Sqrt[1 - L^2/r^2*(1 - 1/r)], r] + c, where c is a constant. The answer is lengthy, but can be plotted. Also, the turning point is given by Solve[1 - L^2/r^2*(1 - 1/r) == 0, r] // First. Apr 13 '17 at 23:23
• Ideally I would want $r$ as a function of $\lambda$ rather than the other way around, is there a simple way of inverting that result? Apr 13 '17 at 23:27
• If you mean in terms of known functions, probably not. So, in the final analysis, the solution as a function of λ must be determined numerically. The advantage of using a symbolic solution at the outset is that numerical integration errors do not accumulate as λ increases. Apr 13 '17 at 23:28
• I just mean in principle, not necessarily in terms of known functions. Interpreting it as a differential equation rather than integrating to find the solution gives me $r$ as a function of $\lambda$, even if it may be a bit clunky. Apr 13 '17 at 23:29

You could differentiate to make a second-order equation without square root problems. Actually, it moves the square trouble to the initial conditions, but that is easier to solve.

foo = D[r'[λ]^2 == 1 - L^2/r[λ]^2*(1 - 1/r[λ]), λ] /. Equal -> Subtract // FactorList
(*  {{-1, 1}, {r[λ], -4}, {r'[λ], 1}, {-3 L^2 + 2 L^2 r[λ] - 2 r[λ]^4 (r^′′)[λ], 1}}  *)

ode = foo[[-1, 1]] == 0     (* pick the right equation by inspection *)
(*  -3 L^2 + 2 L^2 r[λ] - 2 r[λ]^4 (r'')[λ] == 0  *)

icsALL = Solve[{r'[λ]^2 == 1 - L^2/r[λ]^2*(1 - 1/r[λ]), r[0] == 1000} /. λ -> 0,
{r[0], r'[0]}]
(*
{{r[0]  -> 1000,
r'[0] -> -(Sqrt[1000000000 - 999 L^2]/(10000 Sqrt[10]))}, (* negative radical *)
{r[0]  -> 1000,
r'[0] -> Sqrt[1000000000 - 999 L^2]/(10000 Sqrt[10])}}    (* positive radical *)
*)

ics = {r[0], r'[0]} == ({r[0], r'[0]} /. First@icsALL); (* pick negative solution *)

soln = ParametricNDSolve[
{ode, ics},
r, {λ, 0, 2000}, {L}];
Plot[r[30][λ] /. soln, {λ, 0, 2000}]


• Related (differentiating to convert DAE to an ODE): (67642), (91690) Apr 14 '17 at 23:03

As I remarked in earlier comments, this problem can be solved symbolically, although with more difficulty than I had anticipated. The differential equation is given by

r'[λ]^2 == 1 - L^2/r[λ]^2*(1 - 1/r[λ])


which can be rewritten as

Simplify /@ %
(* r'[λ]^2 == (L^2 - L^2 r[λ] + r[λ]^3)/r[λ]^3 *)


It is useful to examine the turning points of this ODE, given by

(L^2 - L^2 r[λ] + r[λ]^3)/r[λ]^3 == 0


Although Solve can handle this equation easily, it is more compact to express the turning points as Root[L^2 - L^2 #1 + #1^3 &, n], where n is 1, 2, or 3. Plotting these roots then gives

Plot[Evaluate@Table[Root[L^2 - L^2 # + #^3 &, n], {n, 3}], {L, 0, 10}, AxesLabel -> {r, L}]


Clearly, the upper (green) turning point curve is desired. It has a lower limit of {L -> (3 Sqrt[3])/2, r -> 3/2}.

Because the ODE itself represents two first order autonomous equations, either can be solved by direct integration.

ss = Integrate[1/Sqrt[1 - L^2/r^2*(1 - 1/r)], r]


I have found no set of Assumptions which causes this expression to yield the desired solution branch. (Edit: Needed Assumptions found.) However, Integrate can treat correctly the more general expression

s = Integrate[r^(3/2)/Sqrt[(r - a) (r - b) (r - c)], r, Assumptions -> a < b < c < r];


although the answer is too long to reproduce here. Given the general answer, the specific answer desired is obtained by

ss = s /. {a -> Root[L^2 - L^2 #1 + #1^3 &, 1], b -> Root[L^2 - L^2 #1 + #1^3 &, 2],
c -> Root[L^2 - L^2 #1 + #1^3 &, 3]};


Its plot for L == 30 (used in the two answers earlier determined numerically) is

With[{L0 = 30}, Quiet@ParametricPlot[{{Chop[-ss /. L -> L0] + 1000, r},
{Chop[ss /. L -> L0] + 1000, r}}, {r, Evaluate@Root[L0^2 - L0^2 # + #^3 &, 3], 1000},
AspectRatio -> 1/GoldenRatio, PlotStyle -> Blue, AxesLabel -> {λ, r}]]


The turning point is 29.4869.

Addendum: L just above its minimum value

The minimum value of L, given above, is

N[3 Sqrt[3]/2]
(* 2.598076211353316 *)


to machine precision. The trajectory for L = 2.5981 is

With[{L0 = 2.5981}, Quiet@ParametricPlot[{{Chop[-ss /. L -> L0], r},
{Chop[ss /. L -> L0], r}}, {r, Evaluate@Root[L0^2 - L0^2 # + #^3 &, 3], 10},
AspectRatio -> 1/GoldenRatio, PlotStyle -> Blue, AxesLabel -> {λ, r},
AxesOrigin -> {0, 0}]]


with a turning point of 1.50372. I presume that the trajectory would skim r == 3/2 for longer and longer times as L further approached 3 Sqrt[3]/2.

Using SolveDelayed -> True inside NDSolve will address your issue. Don't worry about the red color of this addition, it will still work.

f[r_] := 1 - 1/r;

soln = ParametricNDSolve[{r'[λ]^2 == 1 - L^2/r[λ]^2*f[r[λ]], r[0] == 1000}, r,
{λ, 0, 2000}, {L}, SolveDelayed -> True];

Plot[r[30][λ] /. soln, {λ, 0, 2000}, Frame -> True, PlotStyle -> Red]


A new version of SolveDelayed -> True is Method -> {"EquationSimplification" -> "Residual"}.

• Thank you! What does SolveDelayed do? Apr 14 '17 at 7:02
• @qm-arv Check this community.wolfram.com/groups/-/m/t/412961
– zhk
Apr 14 '17 at 7:59
• @qm-arv "EquationSimplification" -> "Residual" solves the ODE as a DAE; it is restricted to the IDA solver, which can handle only machine-precision, real IVPs (no BVPs) at present. Apr 14 '17 at 12:56
• Also, how does Mathematica know too use r'[0] < 0 in this solution? This seems like the "cleanest" solution to the problem if only this were clear to me. Apr 15 '17 at 8:15
• @qm-arv I don't understand what you are asking. Are you looking to find D[r[30][λ], λ] /. soln /. λ -> 0?
– zhk
Apr 15 '17 at 8:19