Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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?

share|improve this question
    
$J_0(x)=cos(xsin(\theta))$? –  Paul Jan 24 '13 at 0:16
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 0:28
    
Ah, I see now... I misread it the first time. –  Paul Jan 24 '13 at 0:34
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
show 1 more comment

1 Answer 1

up vote 3 down vote accepted

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.

share|improve this answer
    
I suspect you're right, my own tests have shown similar results. –  David Z Feb 3 '13 at 22:32
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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