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$.