I was calculating gamma functions in Mathematica while it does not give me an agreed answer.

By definition, $\Gamma[\alpha]=\int_0^\infty t^{\alpha-1}e^{-t}dt$, $\Gamma[\alpha,z]=\int_z^\infty t^{\alpha-1}e^{-t}dt$ and $\Gamma[\alpha,z_1,z_2]=\int_{z_1}^{z_2}t^{\alpha-1}e^{-t}dt$. So we should have $\Gamma[\alpha]-\Gamma[\alpha,z]$ should be equal to $\Gamma[\alpha,0,z]$. But this is not the case. I simply let $\alpha=200$ and $z$ from 1 to 5, he numerical values do not fit; see below.

 α = 200.;
 Table[{Gamma[α] - Gamma[α, z], Gamma[α, 0, z]}, {z, 1, 5}]

{{0.*10^359, 0.00184859}, {0.*10^359, 1.0983*10^57}, {0.*10^359,6.71225*10^91},
 {0.*10^359, 2.41279*10^116}, {0.*10^359,  2.14999*10^135}}

Why is that?


This is precisely the reason Mathematica supplies the three-argument form of Gamma[] (and GammaRegularized[] as well); in computing Gamma[α] - Gamma[α, z] for large α, you are effectively subtracting two large numbers to get a (comparatively) tiny result. (Similar remarks hold for computing Cosh[x] - Sinh[x] just to evaluate Exp[-x]; attempts to perform this for sufficiently large x are doomed.)

Thus, whenever one has to compute differences of incomplete gamma functions (as in this recent example), think about whether it can be reformulated in terms of three-argument Gamma[].

See this related thread on Erf[] as well.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.