I have a function that I need to numerically integrate from $0$ to $+\infty$, given by:

$$I = \int_0^{+\infty} \mathrm{d}x\,x\,T^2(x)f(x)$$

where $T^2$ is an interpolated function that goes to $1$ at small $x$ and decays as $x^{-4}$ at large $x$.

$f$ is an oscillating function with known zeroes:

$$f(x) = j^2_1(\alpha x)\times\cos(\beta\log x)$$

where $j_n$ is the spherical Bessel.

Since I know the zeros $z_n$ of the integrand (thus having control on its oscillations), is it a good idea to split the integral in sub-integrals (from $z_n$ to $z_{n+1}$) and sum all the contributions?

My idea was to change variables ($x = e^u$) and do the sub-integrals with SciPy's scipy.integrate.quad.

  • 2
    $\begingroup$ Hace you checked this post or this post? $\endgroup$
    – nicoguaro
    Jun 19 '17 at 22:52
  • $\begingroup$ I think you could even throw this function into scipy.integrate.quad as it is. The function is smooth and doesn't have singularities (given that $T^2(x)$ is smooth and non-singular), so the quadrature should converge nicely. $\endgroup$ Jun 27 '17 at 7:48
  • $\begingroup$ @Henri Menke: I tried doing that: scipy.integrate.quad raised some "RuntimeWarning: divide by zero encountered" errors, which I assume come from the lower limit of the integral (and the Bessel function). The strange thing is that Mathematica NIntegrate doesn't raise any error and computes the integral over the full range. $\endgroup$
    – HR_8938
    Jun 28 '17 at 8:41

Your Answer

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

Browse other questions tagged or ask your own question.