# Is NDSolve being precise here? Light deflection angle

I am currently trying to calculate deflection angles for the path of photons near Schwarzschild black holes, using the formula

$$\Delta\phi=2\int_{r_0}^\infty \frac{dr}{r^2\sqrt{1/b^2-1/r^2(1-1/r)}}$$

To do this, I solve the equation of motion of the photon numerically:

$$u''(\phi)+u(\phi)=\frac{3}{2}u(\phi)^2$$

(in units where $$2GM=1$$, $$c^2=1$$ and hence $$r_s=1$$ (the Schwarzschild radius).

b = 2;
sol = NDSolveValue[{u''[p] + u[p] == 3/2*u[p]^2, u[0] == 0,
u'[0] == -1/b}, u, {p, 0, 5}]


Then to find $$r_0$$, the radius of closest approach, I use NSolve[sol'[p] == 0, p], which returns {{p -> 1.26921}}. And finally, the integral is evaluated with

2*NIntegrate[ 1/( r^2*Sqrt[1/b^2 - 1/r^2 (1 - 1/r)]), {r,-1/sol[1.2692092946080926], Infinity}]

Which comes out to 1.8498 radians. This seems incorrect, because as shown in the screenshot below from these lecture notes (p.35, author mis-typed his lower integration bound on first integral), the angle should be at least larger than $$\pi$$. Where is the issue?

My only suspicion is the initial condition $$u'(0)=-1/b$$ but it seems correct given the following argument:

$$du/d\phi= (dr/dl)(du/dr)(dl/d\phi)$$

Where $$l$$ is the affine parameter $$dr/dl=1$$, the speed of light (this is where I'm not too convinced as affine parameter is not proper time) then use E-L equations to get $$dl/d\phi$$ and $$u=1/r$$.

• Equation is not needed in this problem. The angle is determined by the integral. – Alex Trounev Dec 4 '18 at 0:36
• @AlexTrounev, completely agree that angle is determined by the integral, the equation (from lecture notes) I just included so that we can verify the answer Mathematica gives. My point is that the equation implies that the angle (given by the integral) must be greater than pi, but Mathematica gives a number smaller than pi. – tr416 Dec 5 '18 at 8:12
• This should be exactly solvable in terms of Weierstrass elliptic functions; if you are interested in that, I can try to write something up. – J. M. is computer-less Jan 5 at 8:00

The angle is calculated as an integral depending on the parameter b. Calculate the minimum distance

 r0 = r /. First[Solve[1/b^2 - 1/r^2 (1 - 1/r) == 0, r]]

(*Out[]= ((2/3)^(1/3) b^2)/(-9 b^2 + Sqrt[3] Sqrt[27 b^4 - 4 b^6])^(
1/3) + (-9 b^2 + Sqrt[3] Sqrt[27 b^4 - 4 b^6])^(1/3)/(
2^(1/3) 3^(2/3))*)


Calculate the angle for example when b = 4

 2*NIntegrate[
1/(r^2*Sqrt[1/b^2 - 1/r^2 (1 - 1/r)] /. b -> 4.), {r,
Re[r0 /. b -> 4.], Infinity}]


(* Out[]= 4.00032*)

In this problem there is a minimum value of b=b0=(3 Sqrt[3])/2. For b < b0 the minimum radius becomes complex, for example for b = 2

r0 /. b -> 2.
(*1.191487883953119 - 0.5088517788327379 I*)


This means that light is captured by a black hole.