1
$\begingroup$

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 Edit: adding equation source

$\endgroup$
4
  • $\begingroup$ Should $x$ be $u$ in that arctan function? $\endgroup$
    – Charlie S
    Oct 9 '20 at 18:51
  • $\begingroup$ Please make sure the question is correctly formulated. Do you integrate over $u$ or over $x$? $\endgroup$ Oct 10 '20 at 5:45
  • $\begingroup$ Thank you. After I read the cited reference (Smith, 1937), I have found there are some errors in the above pasted source paper, and I have assured that the intergration is over x. I have corrected the above equation and code. However, the same warning is there. $\endgroup$
    – NoVel
    Oct 10 '20 at 14:09
  • $\begingroup$ Have you tried to split the integral at the sum? I feel like the first summand could have an analytical solution, simplyfing the task and numerical effort $\endgroup$
    – MPIchael
    Oct 15 '20 at 14:35
1
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ in addition the first integral has an analytical solution, so they could remove this first. $\endgroup$
    – porphyrin
    Oct 14 '20 at 8:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.