Upper incomplete Gamma function numerical differences

Question

I have the following issue. Using Mathematica 11.2 (with Rubi loaded) I find that e.g.
In[504]:= Gamma[-1/3, 2] // N
Out[504]= 0.0353296 + 5.01821*10^-16 I
The imaginary part is absurd, and the correct approximation is the real part of the output. Indeed, if I look for the result using the Wolfram query, it is the correct one. I am pretty sure that this is a numerical issue since I am getting things like
In[507]:= Gamma[-10/3, 2] // N
Out[507]= 0.00234763 + 0. I
What can I do to resolve this? Is it something known?

CA Trevillian · Accepted Answer

I ran these on the WolframCloud and seemingly the issue is fixed now on v12.2.0.
$Version
Gamma[-1/3,2]
N@%
Gamma[-1`20/3`20, 2`20]

12.2.0 for Linux x86 (64-bit) (November 16, 2020)
Gamma[-1/3,2]
0.0353296
0.035329560217661993

If none of these methods work, I would either do these calculations on the cloud or upgrade your license to the newest version. The issue is one of numerical precision & underlying algorithms.
I will say that the relative magnitude of the complex term is small enough in comparison to the real term to trust in the use of Chop as recommended by others.

Bob Hanlon · Answer

$Version

(* "12.0.0 for Mac OS X x86 (64-bit) (April 7, 2019)" *)

Gamma[-10/3, 2] // N

(* 0.00234763 + 0. I *)

A transformation for Gamma[a, z] is
repl = Gamma[a_, z_] -> (-z^a + E^z*Gamma[1 + a, z])/(E^z*a);

Verifying,
Gamma[a, z] == (Gamma[a, z] /. repl) // FullSimplify

(* True *)

Use ReplaceRepeated
Gamma[-10/3, 2] //. 
   Gamma[a_?Negative, z_] -> (-z^a + E^z*Gamma[1 + a, z])/(E^z*a) // 
  Simplify // N

(* 0.00234763 *)

Upper incomplete Gamma function numerical differences

2 Answers

Add your own answers!

Ask a Question