How can I integrate a function related to the modified bessel function of the first kind?
조회 수: 11 (최근 30일)
이전 댓글 표시
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);
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);
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
채택된 답변
Bruno Luong
2022년 9월 6일
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)
참고 항목
카테고리
Help Center 및 File Exchange에서 Bessel functions에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!