problem with gammainc function

neo
2012 12 18
2 个回答
回答已采纳
4 次查看（30 天）

0 个投票

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...

2 个评论
显示 无 隐藏 无

neo
neo 2012-12-25
编辑：neo 2012-12-25
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

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Walter Roberson
Walter Roberson 2012-12-19

1 个投票

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)

更多回答（1 个）

Peter Perkins
Peter Perkins 2012-12-19

1 个投票

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.

2 个评论
显示 无 隐藏 无

neo
neo 2012-12-25
编辑：neo 2012-12-25
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
Peter Perkins 2012-12-26
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.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by