How can I integrate a function related to the modified bessel function of the first kind?

Question

祥宇崔 2022-9-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1796290-how-can-i-integrate-a-function-related-to-the-modified-bessel-function-of-the-first-kind

评论：祥宇崔 2022-9-6

采纳的回答： Bruno Luong

在 MATLAB Online 中打开

Thank you for reading my question!

The thing is that I'm tring to integrate a rice PDF, and it should be 1.

beta = 5 ;
sigma = 1.4/3;
distribution_func = @(x) x./sigma.^2.*exp(-(x.^2+beta.^2)./2./sigma.^2).*besseli(0,beta.*x./sigma.^2);
integral(distribution_func,0,100);
Warning: Inf or NaN value encountered.

As you can see, due to

can go to Inf in a short range, besseli function can't be used.

And I think

might be big alone, but it would be small with the scaling in distribution_func. So I try to constuct my own bessel function and it's like:

beta = 5 ;
sigma = 1.4/3;
integ_func = @(theta,x) x./sigma.^2./pi.*exp(-(x.^2+beta^2)/2/sigma^2+beta*x/sigma^2.*cos(theta));
distribution_func = @(x) integral(@(theta) integ_func(theta,x),0,pi);
integral(distribution_func,0,1);
Error using integralCalc/finalInputChecks
Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the 'ArrayValued' option to true.

Error in integralCalc/iterateScalarValued (line 315)
                finalInputChecks(x,fx);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 87)
Q = integralCalc(fun,a,b,opstruct);

The equation is right actually, but it seems that Integral can't be used in this way.

So I try trapz:

beta = 5 ;
sigma = 1.4/3;
integ_func = @(theta,x) x./sigma.^2./pi.*exp(-(x.^2+beta^2)/2/sigma^2+beta*x/sigma^2.*cos(theta));
distribution_func = @(x) integral(@(theta) integ_func(theta,x),0,pi);
x = 0:0.01:1e3;
y = zeros(1,length(x));
for i = 1:length(x)
    y(i) = distribution_func(x(i));
end
disp(trapz(x,y));

It' ok, but I don't like it, I think it's clumsy and inaccurate under certain situations.

The last method I may try is double integral since there is two integration in my equation.

Could you give me some advice about some easy and accurate ways to reach my goal?

Thanks in advance!!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Bruno Luong 2022-9-6

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1796290-how-can-i-integrate-a-function-related-to-the-modified-bessel-function-of-the-first-kind#answer_1043535

在 MATLAB Online 中打开

You can also "vectorize" distribution_func by applying arrayfun

beta = 5 ;
sigma = 1.4/3;
integ_func = @(theta,x) x./sigma.^2./pi.*exp(-(x.^2+beta^2)/2/sigma^2+beta*x/sigma^2.*cos(theta));
distribution_func = @(x) arrayfun(@(x)integral(@(theta) integ_func(theta,x),0,pi), x);
integral(distribution_func,0,1)

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

祥宇崔 2022-9-6

Thanks!

I'll try yo figure out what is Arrayfun.

请先登录，再进行评论。

Answer 2

祥宇崔 2022-9-6

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1796290-how-can-i-integrate-a-function-related-to-the-modified-bessel-function-of-the-first-kind#answer_1043510

在 MATLAB Online 中打开

It seems that this works:

beta = 5 ;
sigma = 1.4/3;
integ_func = @(theta,x) x./sigma.^2./pi.*exp(-(x.^2+beta^2)/2/sigma^2+beta*x/sigma^2.*cos(theta));
distribution_func = @(x) integral(@(theta) integ_func(theta,x),0,pi,'ArrayValued',true);
integral(distribution_func,0,10);