15

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?

J. M.
  • 3,155
  • 28
  • 37
David Z
  • 3,383
  • 2
  • 27
  • 34
  • 1
    Integrated over $0\le\theta \le\pi$... It's one of the integral definitions of the Bessel function. – David Z Jan 24 '13 at 00:28
  • 4
    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)]]$. – Stefano M Jan 24 '13 at 16:21
  • @StefanoM yes, that's right. – David Z Jan 24 '13 at 19:57
  • 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. – J. M. May 14 '13 at 04:30

1 Answers1

3

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.

Ondřej Čertík
  • 2,930
  • 18
  • 40