Bessel function solution: Mathmatical point of view
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
0 commentaires
Réponses (1)
Daniel Baboiu
le 3 Nov 2011
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
3 commentaires
Voir également
Catégories
En savoir plus sur Bessel functions dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!