Bessel function solution: Mathmatical point of view

Question

Amin 2011-11-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/20247-bessel-function-solution-mathmatical-point-of-view

在 MATLAB Online 中打开

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Daniel Baboiu 2011-11-3

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/20247-bessel-function-solution-mathmatical-point-of-view#answer_26761

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