How can I integrate this function?

Hey guys, I have this function here: where is the modified Bessel function of the first kind of order K-1. I want to integrate from b (which is a constant I decide on) to infinity. And K and B are also constants which I decide on. This integral that I want to perform is insanely close to the Marcum-Q integral which MATLAB already has a function for (https://www.mathworks.com/help/signal/ref/marcumq.html). Unfortuantely though, I don't think its possible to use it because my is very slightly different. I also tried developing my own function and failed. Any help is appreciated.

Answer 1

[Star Strider]

The function as presented appears to be a bit ambiguous, so I’m not certain which (if either) of these is correct.

That aside, the result is NaN with an infinite upper limit for both of them —

p = @(r,K,B) 0.5*(r./(2*K*B)).^((K-1)/2) .* exp(-(r+2*K*B)/2) .* besseli(K-1,sqrt(2*K*B*r));

K = 2;
B = 1;
b = rand

b = 0.0492

Result = integral(@(r)p(r,K,B), b, 1E5)

Result = 1.0000

.

[Torsten]

@Star Strider

The first functional form is correct when comparing with the Marcum-Q function. But the 1/2 must be left outside the exponentiation and the r inside.

p = @(r,K,B) (r./(2*K*B)).^((K-1)/2) .* exp(-(r+2*K*B)/2) .* besseli(K-1,sqrt(2*K*B*r));
K = 2;
B = 1;
b = rand

b = 0.0449

b = 0.6997

b = 0.6997

Result = 0.5*integral(@(r)p(r,K,B), b, 1E5)

Result = 0.9918

[Walter Roberson]

besseli is non-decreasing for odd n, and when K=2 then n=K-1 happens to be odd. Under that circumstance you can show that the function value goes to a constant as r goes to infinity, and infinite integral of a nonzero constant is infinite.

If K happens to be odd then there could be a different outcome.

[Star Strider]

@Torsten — Thank you for the clarification. Code changed.

It still gives a NaN result with an infinte upper integration limit, so either I must have not chosen the parameter values correctly (I have no idea what they should be, however I wanted to test the code), or the integral is simply undefined, regardless.

[Ali Almakhmari]

Thanks guys. This gives me a solid ground to do what I wanted to originally do.

[Star Strider]

As always, our pleasure!

[Torsten]

@Star Strider

Note that the r is still outside the exponentiation in your definition of p.

[Star Strider]

Fixed. I was staring so intently at the LaTeX expression that I overlooked that.

Answer 2

[Walter Roberson]

Assuming the the posted formula is correct in having the exp() raised as an exponent to the (r/2KB) then:

Because you want to integrate to r = infinity, examine your p(r) as r increases towards infinity.

exp(-(r+2KB)/2) as r increases towards infinity heads to exp(-infinity) which head towards 0 from above.

That 0-limit is being multiplied by a constant (K-1)/2 which continues to give the 0-limit.

That 0-limit is the exponent of r/(2*K*B) and so that limit would be r^0 which would be limit-1 -- a non-zero constant

The leading 1/2 is a positive multiplicative constant, and a positive constant times a limit-constant gives a positive constant in the limit.

Thus as r goes to infinity the part of the expression before the besseli goes to positive constant.

Now, the integration of a positive constant out to infinity is infinity, so if we ignore the besseli for the moment, we see we are at risk of an infinite outcome of the integral.

When would the true integral potentially be below infinity? Only if the besseli() portion goes towards zero. But for real-valued constant nu, besseli(nu, x) increases without bound with increasing x. And sqrt(2*K*B*r) is strictly increasing in increasing r.

Therefore, if the exp() is properly part of the exponent of (r/2KB) then the integral of p as r goes to infinity is going to be infinite.

... Now if the exp() should not be part of the exponent, if it should be at the same baseline as the (r/2KB) then the analysis would be notably different.

[Torsten]

Assuming the the posted formula is correct in having the exp() raised as an exponent to the (r/2KB)

It doesn't seem to be the case:

https://www.mathworks.com/help/signal/ref/marcumq.html