How to numerically integrate bessely function

1 次查看（过去 30 天）

显示更早的评论

David Zhang 2012-2-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/29920-how-to-numerically-integrate-bessely-function

编辑： Cedric 2013-10-20

在 MATLAB Online 中打开

Hi All, i want to get the integral result for bessely function as

int{a}{b}int{a}{b}bessely(1, x1-x2)dx1dx2

but when numerically calculate that using Gaussian quadrature points, the values in bessely(1, 0) is a infinite number, and the limit as limit{x->0}bessely(1, x) are not existed. But when i used the quad2d to calculate this integral, a correct answer can be obtained. So i like to know how the results can be numerically obtained? Thanks.

My code is as

clear
clc
close all
a=-8;
b = 6;
% covpara=1;
%%number of quadrature point
n=50;
% % % % the whole gaussian quadrature point
[Q, W]=LegGaus(n,a,b);
covfun_A=@(x1,x2) abs(x1-x2).*bessely(1,abs(x1-x2));
Matrix_A = bsxfun(covfun_A,Q.',Q); % f = @(x,y)....
%A_1 impose the diagnoal elements as a very small number, i.e., bessely 
Matrix_A(1:length(Matrix_A)+1:length(Matrix_A)*length(Matrix_A))=0; 
%A_1 is the numerical result using Gaussian quadrature, using tensor
%product of 1D to get 2D result
A_1 =(W)'*Matrix_A*(W);
%A_2 is the 'exact' result from matlab
A_2= (quad2d(@(x1,x2) (abs(x1-x2).*bessely(1,abs(x1-x2))),a,b,@(x1)x1,b) +...
    quad2d(@(x1,x2) (abs(x1-x2).*bessely(1,abs(x1-x2))),a,b,a,@(x1)x1));
A_2/A_1