Numerical integration problem: IntegrationWarning The integral is probably divergent, or slowly convergent

Question

I am trying to get the numerical integration of a function using scipy's integrate.quad as follows.

$$ \begin{equation} G (\alpha) = \frac{4\alpha}{\pi}\int_0^{\infty} x e^{-\alpha x^2} {(\pi/2}+\mathrm{tan}^{-1}[Y_0(x)/J_0 (x)]) \mathrm{d}x \end{equation} $$

import numpy as np
from scipy import integrate
from scipy.special import k0,j0,y0,k1
def G(alpha=743711.5,T=5.5e-5,sw=10.65):
    pi = 3.14
    fun = lambda x: x*np.exp(-1*alpha*x**2)*(pi/2+np.arctan(y0(x)/j0(x)))
    val,err = integrate.quad(fun,0, np.inf)
    return val,err
val,err = G()
print (val,err)

However, I get "The integral is probably divergent, or slowly convergent." I have tried to set a very large limit, such as limit = 10000000. However, the same warning is. Does anyone know how to solve the problem?

https://doi.org/10.1029/TR033i004p00559

Answer 1

The problem is probably that the value of $\alpha$ is large compared to the values of $x^2$ that the numerical integrator is probably choosing. For any value $x \gg 1/\sqrt{\alpha}$ the integrand will be close to zero (for large enough values it will be exactly zero in floating point). The integrand is also zero at $x = 0$. So there is no point in setting the limit to a large value. On the contrary, try setting it to something like $10/\sqrt{\alpha}, 20/\sqrt{\alpha}, ...$ and check for convergence of the results.

Perhaps best would be to define a new variable $u = x/\sqrt{\alpha}$ and integrate with respect to $u$. Then the numerical integrator should probably work.