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

16 views (last 30 days)
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!!

Accepted Answer

Bruno Luong
Bruno Luong on 6 Sep 2022
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)

More Answers (2)

祥宇 崔
祥宇 崔 on 6 Sep 2022
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);

Sam Chak
Sam Chak on 6 Sep 2022
Have you checked if this besseli() function for Modified Bessel function of first kind is applicable for your case?

Products


Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!