problem with gammainc function

I am calculating incomplete gamma function using gammainc gammainc(0.04344,0,'upper') and I get the result as zero

but when I calculate same integral numerically I get result 2.6. I have also verified this results with other software, which gives the same result.

for example (using MuPAD)

double(feval(symengine,'igamma', 0, 0.043444))

gives me the result 2.6

Why I am not able to get same result with gammainc...

[Jonathan Sullivan]

Did you?

Both MATLAB and WolframAlpha agree.

http://www.wolframalpha.com/input/?i=IncompleteGamma%5B0.04344%2C0%5D

[neo]

Thanks Jonathan, but actually they don't agree.

Try Gamma[0,0.43444] in Mathematica it gives results 2.6 In my problem a=0 (and x is any positive value)

Try this link http://www.wolframalpha.com/input/?i=Gamma%5B0%2C0.04344%5D

Answer 1

[Walter Roberson]

If you go by the integrals on the gammainc() documentation page and examine them with respect to approaching a = 0, you find that you have a term which is 1 / gamma(0); and since gamma(0) is infinity, 1 / gamma(0) approaches 0, which could lead you to expect that gammainc(x,0) (using lower tail) should be 0.

But when you continue the analysis, you realize that that is only true provided that the integral itself does not approach infinity.... which it does. Thus for the case a = 0, the function must proceed either by way of deeper analysis, or by way of approximation.

The symbolic math igamma() function proceeds by way of deeper analysis and converts igamma(0,x) to Ei(x), the exponential integral.

Do we know what gammainc() does? Not for sure. But an important clue is the statement in the documentation that,

For small X and A, gammainc(X,A) is approximately equal to X^A, so gammainc(0,0) = 1.

If we substitute in A=0 for that, we see that gammainc() plausibly considers gammainc(X,0) to be approximately equal to X^0 which is 1. As that is the lower tail, the upper tail would then be approximated as 1 minus that 1, giving 0. And that fits the actual returned value for gammainc(X,0,'upper'), so it probably isn't worth continuing with guessing.

So... in this boundary case, gammainc() is likely approximating and getting the wrong approximation, whereas igamma() knows the correct analytic conversion.

If you know ahead of time that your A will be 0, you should probably use the exponential integral, expint(x)

Answer 2

[Peter Perkins]

This is just a standardization issue. The doc for MATLAB's gammainc and MuPAD's igamma clearly show that the former is standardized by 1/gamma(a), while the latter is not. And so (e.g. for a == .1):

>> gammainc(0.043444,.1,'upper') * gamma(.1)
ans =
          2.23413826678516
>> double(feval(symengine,'igamma', .1, 0.043444))
ans =
          2.23413826678516

You can't do that comparison at a==0, but in the limit, gamma(0) is Inf and expint(0.043444) is finite, so gammainc(0.043444,0,'upper') ought to be 0.

[neo]

Thanks peter, but I was looking for the case where a=0

So I think I should use expint(x) function as suggested by Walter or use MUPAD. I also got results by this code also

E1=@(z)exp(-z)./z;

quadgk(E1,0.043444,Inf)

[Peter Perkins]

neo, you don't seem to have gotten my point. Read the descriptions of the functions you are calling, and you will find they are not calculating the same things. There is not "a problem" with MATLAB's gamminc when a is 0. It may be that what you want is in fact expint(0.043444), but gammainc(0.043444,0,'upper') is correct when it returns 0.