CODEWITHSUNDEEP
matlabintegralbessel-functions
Nikolai Liubimov
Nikolai Liubimov

Reputation: 76

Integration of 0-th order Bessel function using MATLAB

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?

Upvotes: 3

Views: 4205

Answers (2)

Rody Oldenhuis
Rody Oldenhuis

Reputation: 38042

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.

Upvotes: 1

Dan
Dan

Reputation: 45752

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)

Upvotes: 1

Related Questions