Finding Scattering Angle by solving Integral numerically with Singularities

Question

I was trying to do a computational project from Computational Physics by Steven Koonin. The problem requires us to solve two integrals each of which has singularity at the lower limit and numerically i have been not able to solve this problem. The integral is $$ \theta = 2b \left [ \int_b^{r_{max}} \frac{dr}{r^2}(1- \frac{b^2}{r^2})^{-1/2} - \int_{r_{min}}^{r_{max}} \frac{dr}{r^2}(1-\frac{b^2}{r^2} - \frac{V}{E})^{-1/2} \right ] $$

Here, $r_{min}$ is the root of argument in second square root and b is root for first integral.

Can anyone suggest some method to solve this numerically preferably by simpsons 1/3rd or trapezoidal rule? I am working with python.

Answer 1

You cannot use Simpson or Trapezoidal rule for solving these integrals because both methods utilize the two bounds when evaluating the functions. In fact, most quadrature methods you might try first will do the same. Hence, you will get NaN or "div by zero" errors.

What you need instead is a quadrature that does not utilize the bounds for integration, such as Gaussian quadrature. Or, if you have access to it, numerical schemes that omit problematic points (see NAG's d01slc for instance).

Note, though, that if you let $\mu=b/r$, then you end up with an integral like, $$ I_1\sim\int \mathrm d\mu\frac{1}{\sqrt{1-\mu^2}} $$ which looks a lot like a Chebyshev-Gauss quadrature of the first kind, which might be the direction to investigate first. Transforming the second integral likely will be a bit uglier, but can be done as well.