15
$\begingroup$

I've heard anecdotally that when one is trying to numerically do an integral of the form

$$\int_0^\infty f(x) J_0(x)\,\mathrm{d}x$$

with $f(x)$ smooth and well-behaved (e.g. not itself highly oscillatory, nonsingular, etc.), then it will help accuracy to rewrite it as

$$\frac{1}{\pi}\int_0^\pi \int_0^\infty f(x) \cos(x\sin\theta) \,\mathrm{d}x\,\mathrm{d}\theta$$

and perform the inner integral numerically first. I can't see any reason I should expect this to work, but then again the accuracy of a numerical method is rarely obvious.

Of course I know the best way to actually do it is to use a method optimized for oscillatory integrals like this, but for curiosity's sake, suppose I limit myself to using some quadrature rule. Can anyone confirm or refute that making this transformation tends to improve the accuracy of the integral? And/or point me to a source explaining it?

$\endgroup$
4
  • 1
    $\begingroup$ Integrated over $0\le\theta \le\pi$... It's one of the integral definitions of the Bessel function. $\endgroup$
    – David Z
    Jan 24, 2013 at 0:28
  • 4
    $\begingroup$ So your question is: Given generic $N$-point quadrature formulas $Q_\infty^N[\cdot]$ on $[0,\infty)$ and $Q_\pi^N[\cdot]$ on $[0,\pi]$, is $Q^{N\cdot M }_\infty[f\,J_0]$ worse or better than $Q^M_\pi[Q_\infty^N[f(x)\,\cos(x\sin\theta)]]$. $\endgroup$
    – Stefano M
    Jan 24, 2013 at 16:21
  • $\begingroup$ @StefanoM yes, that's right. $\endgroup$
    – David Z
    Jan 24, 2013 at 19:57
  • $\begingroup$ FWIW, one of the most efficient methods for evaluating the zeroth-order Bessel function is the trapezoidal rule, which is well-known to give very accurate results when integrating periodic integrands over one period (even better than the usual standard, Gaussian quadrature). So: it might help, it might not. $\endgroup$
    – J. M.
    May 14, 2013 at 4:30

1 Answer 1

3
$\begingroup$

I don't think it makes any difference. You have to choose high enough quadrature for the integral over $\theta$ so that it is equal to the Bessel function $J_0$. I chose order 20 in the example below, but you always have to do convergence with regards to the exact function and interval that you integrate over. Then I did convergence with $n$, the order of the Gaussian quadrature of the integral over $x$. I chose $f(x) = e^{-x} x^2$ and use domain $[0, x_\text{max}]$, you can change $x_\text{max}$ below. I got:

 n      direct         rewritten
 1  0.770878284949  0.770878284949
 2  0.304480978430  0.304480978430
 3  0.356922151260  0.356922151260
 4  0.362576361509  0.362576361509
 5  0.362316789057  0.362316789057
 6  0.362314010897  0.362314010897
 7  0.362314071949  0.362314071949
 8  0.362314072182  0.362314072182
 9  0.362314072179  0.362314072179
10  0.362314072179  0.362314072179

As you can see, for $n=9$ both integrals are fully converged to 12 significant digits.

Here is the code:

from scipy.integrate import fixed_quad
from scipy.special import jn
from numpy import exp, pi, sin, cos, array

def gauss(f, a, b, n):
    """Gauss quadrature"""
    return fixed_quad(f, a, b, n=n)[0]

def f(x):
    """Function f(x) to integrate"""
    return exp(-x) * x**2

xmax = 3.

print " n      direct         rewritten"
for n in range(1, 20):
    def inner(theta_array):
        return array([gauss(lambda x: f(x) * cos(x*sin(theta)), 0, xmax, n)
            for theta in theta_array])
    direct = gauss(lambda x: f(x) * jn(0, x), 0, xmax, n)
    rewritten = gauss(inner, 0, pi, 20) / pi
    print "%2d  %.12f  %.12f" % (n, direct, rewritten)

You can play with this yourself, just change xmax, possibly you might need to split the interval $[0, \infty]$ into elements and integrate element by element. You can also change the function f(x). Make sure that you always converge the integral rewritten = gauss(inner, 0, pi, 20) / pi, i.e. start with some low order and keep increasing it until the printed results stop changing.

$\endgroup$
1
  • $\begingroup$ I suspect you're right, my own tests have shown similar results. $\endgroup$
    – David Z
    Feb 3, 2013 at 22:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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