Bessel function solution: Mathmatical point of view

Hi,

I am integrating a bessel function, two ways. But the results are different.

1) using quad,like this:

q1 = quadl( @(u)testf_bessel(u),0, 100);

where 'testf_bessel) is a function:

function arg = testf_bessel(x)
arg = double((bessel(1.5,x)).^2./x.^3);

2) first

besselj(1.5,x).^2 ./ x.^3

then I get : (2*(cos(x) - sin(x)/x)^2)/(pi*x^4) now integrate:

q2 = quad(@(x) ( (2*(cos(x) - sin(x)./x).^2)./(pi*x.^4) ),0,100);

in the range of 0 to 100 they give the same result. But if I use a for loop and change the high-limit of integral from 10 to 10000,and plot it, I see that q1 is always at 0.133, but q2 drops to zero after some time.

I am not a mathematician, anyone could give me some explanation?

Thanks

Answer 1

[Daniel Baboiu]

This is a problem of the quad and quadl functions.

This phenomenon appears regardless of the use of quad or quadl (there are subtle differences in the adaptive algorithm). Mathematically, besselj(1.5,x) and (2*(cos(x) - sin(x)./x).^2)./(pi*x.^4)are identical, but the implementation is different. This leads to a significant difference in computing the position of the sample points. You have highly oscillating functions over intervals much larger than the period of the oscillation, even though the amplitude decreases quickly; numerical integration algorithms may be unstable in such situations.

If you use the form [q,n]=quadl(...) then n is the number of function evaluations. You also have an additional parameter to quad/quadl to control tolerance:

>>[q,count]=quad(@(x) double( (2*(cos(x) - sin(x)./x).^2)./(pi*x.^4) ),0,100)

q = 0.1333 count = 94

>>[q,count]=quad(@(x) double( (2*(cos(x) - sin(x)./x).^2)./(pi*x.^4) ),0,1000)

q = 4.3969e-07 count = 14

>>[q,count]=quad(@(x) double( (2*(cos(x) - sin(x)./x).^2)./(pi*x.^4) ),0,1000,1e-10)

q = 0.1333 count = 722

>>[q,count]=quad(@(x) double( (2*(cos(x) - sin(x)./x).^2)./(pi*x.^4) ),0,100000,1e-10)

q = 2.2069e-12 count = 14

>>[q,count]=quad(@(x) double( (2*(cos(x) - sin(x)./x).^2)./(pi*x.^4) ),0,100000,1e-20)

Warning: Maximum function count exceeded; singularity likely.

q = 0.1333 count = 10086

[Amin]

I tried with quadgk. I don't know what exactly quadgk is doing, but this was the result!

>>[q,count]=quadgk(@(x) double( (2*(cos(x) - sin(x)./x).^2)./(pi*x.^4) ),0,inf)

q = 0.1333 count = 4.5278e-08

Does this mean I can count on quadgk in my calculations? :-?

thanks for the help, Daniel Baboiu

[Amin]

what does 'inf' mean here?

[Walter Roberson]

inf means infinity in MATLAB.