Integration of 0-th order Bessel function using MATLAB

Question

I have simple question. I'm trying to evaluate improper integral of 0th order Bessel function using Matlab R2012a:

v = integral(@(x)(besselj(0, x), 0, Inf)

which gives me v = 3.7573e+09. However this should be v = 1 in theory. When I'm trying to do

v = integral(@(l)besselj(0,l), 0, 1000)

it results to v = 1.0047. Could you briefly explain me, what is going wrong with integration? And how to properly integrate Bessel-type functions?

Rody Oldenhuis · Answer 1 · 2013-04-05T15:48:42.420

At first I was sceptical that taking an integral over a Bessel function would produce finite results. Mathematica/Wofram Alpha however showed that the result is finite, but it is not for the faint of heart.

However, then I was pointed to this site where it is explained how to do it properly, and that the value of the integral should be 1.

I experimented a bit to verify the correctness of their statements:

    F = @(z) arrayfun(@(y) quadgk(@(x)besselj(0,x), 0, y), z);

    z = 10:100:1e4;

    plot(z, F(z))

which gave:

bessel integral

so clearly, the integral indeed seems to converges to 1. Shame on Wolfram Alpha!

(Note that this is kind of a misleading plot; try to do it with z = 10:1e4; and you'll see why. But oh well, the principle is the same anyway).

This figure also shows precicely what the problem is you're experiencing in Matlab; the value of the integral is like a damped oscillation around 1 for increasing x. Problem is, the dampening is very weak -- as you can see, my z needed to go all the way to 10,000 just to produce this plot, whereas the oscillation amplitude was only decreased by ~0.5.

When you try to do the improper integral by messing around with the 'MaxIntervalCount' setting, you get this:

    >> quadgk(@(x)besselj(0,x), 0, inf, 'maxintervalcount', 1e4) 

    Warning: Reached the limit on the maximum number of intervals in use.
    Approximate bound on error is 1.2e+009. The integral may not exist, or
    it may be difficult to approximate numerically.
    Increase MaxIntervalCount to 10396 to enable QUADGK to continue for 
    another iteration. 
    > In quadgk>vadapt at 317
      In quadgk at 216

It doesn't matter how high you set the MaxIntervalCount; you'll keep running into this error. Similar things also happen when using quad, quadl, or similar (these underly the R2012 integral function).

As this warning and plot show, the integral is just not suited to approximate accurately by any quadrature method implemented in standard MATLAB (at least, that I know of).

I believe the proper analytical derivation, as done on the physics forum, is really the only way to get to the result without having to resort to specialized quadrature methods.

http://www.physicsforums.com/showthread.php?t=570379 looks like it should be 1 :/ see post #9 — Dan, Apr 05 '13 at 14:07
@Dan: [Check mathematicas answer](http://integrals.wolfram.com/index.jsp?expr=BesselJ%5B0%2C+x%5D&random=false); that sure don't look like '1' to me :) (or well, a '2' than, since it's improper, but...same thing) — Rody Oldenhuis, Apr 05 '13 at 14:48
Hmmm, didn't load nicely on my side. This has integral of 1: http://www.wolframalpha.com/input/?i=integrate+first+bessel+function+from+0+to+infinity but it looks like a reflection of yours. So maybe OP is integrating the wrong bessel function? but it doesn't look like it based on the result where `inf` is replaced with `1000`. I wonder if it starts to converge if that `1000` is replaced with `1000000`. Don't have access to matlab now to check though. — Dan, Apr 05 '13 at 15:02
@Dan: Well, see my edit; I was indeed wrong. Interesting problem! I'll experiment a bit more. — Rody Oldenhuis, Apr 05 '13 at 15:12
@Dan: Anyway, the question was about *zeroth* order bessel funtion, so, [use the corresponding mathematica solution](http://www.wolframalpha.com/input/?i=integrate+zeroth+bessel+function+from+0+to+infinity). — Rody Oldenhuis, Apr 05 '13 at 15:13
Thank's! I've realized that it would be better to switch to some analytical expression that could be approximated more efficiently — Nikolai Liubimov, Apr 05 '13 at 16:24

score 1 · Answer 2 · answered Apr 05 '13 at 14:04

From the docs to do an improper integral on an oscillatory function:

q = integral(fun,0,Inf,'RelTol',1e-8,'AbsTol',1e-13)

in the docs the example is

fun = @(x)x.^5.*exp(-x).*sin(x);

but I guess in your case try:

q = integral(@(x)(besselj(0, x),0,Inf,'RelTol',1e-8,'AbsTol',1e-13)

Integration of 0-th order Bessel function using MATLAB

2 Answers2