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

14 views (last 30 days)

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

[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)

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?

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)

Thanks!
I'll try yo figure out what is Arrayfun.

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);

hh, it works!

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?

Hi,Chak!
Maybe not
But I think I can use the scale term in besseli to make it applicable

### Categories

Find more on Bessel functions in Help Center and File Exchange

R2020a

### Community Treasure Hunt

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

Start Hunting!