Bessel function solution: Mathmatical point of view
조회 수: 1 (최근 30일)
이전 댓글 표시
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
답변 (1개)
Daniel Baboiu
2011년 11월 3일
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
참고 항목
카테고리
Help Center 및 File Exchange에서 Bessel functions에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!