Warning when confluent hypergeometric function HypergeometricU is wrong?

Question

From the definition as an integral (HypergeometricU/details), this function must be positive. However, it gives negative numbers in some cases with no warning of a potential error. For example,

In[228]:= HypergeometricU[2., -97., 177.]
Out[228]= -1.18092*10^64

The negative output is not small, so it is not an apparent underflow problem. I can imagine reasons why Mathematica's algorithm gets it wrong, and this error is easy to catch. But I am concerned that it is making errors in other cases. Is there a way to get a warning when the result is not trustworthy?

Answer 1

$Version

(* "12.2.0 for Mac OS X x86 (64-bit) (December 12, 2020)" *)

Clear["Global`*"]


HypergeometricU[2., -97., 177.]

(* -1.18092*10^64 *)

Any calculation done with machine precision is done with the understanding that "you get what you get." Machine precision is fast but neither tracks nor attempts to control precision. Use arbitrary-precision:

HypergeometricU @@ SetPrecision[{2., -97., 177.}, 15]

(* 0.00001308368035127 *)

Or for an exact solution

val = HypergeometricU @@ Rationalize[{2., -97., 177.}]

(* -(1612891645768307438642571599604930828030817964923418757903738393704630414517\
842865730334134639989235034374108215803320414126362646541365911522932125411703\
13058197149329186095569904019297562571667/
    48498254399355628616222215525374070469654954510969291222415885140849998225\
60734976299979989248181888803143680000000000000000000000) + \
(12759366726360372317225946603232298388725889904857641932911677232166739818527\
506078285603785264758197470099770656607352615543190362442650920745312680944219\
211068340938900712983219818678309780201593 E^177 Gamma[0, 177])/
  2155477973304694605165431801127736465317997978265301832107372672926666587804\
771100577768884110303061690286080000000000000000000000 *)

Block[{$MaxExtraPrecision = 100}, N[val, 15]]

(* 0.0000130836803512679 *)

Answer 2

In this case it really helps to use an explicit formula for $U(2,b,z)$, which can be expressed in terms of the much more stable exponential integral function:

U2[b_, z_] = HypergeometricU[2, b, z] // FullSimplify
(*    (1 - E^z (2 - b + z) ExpIntegralE[2 - b, z])/(-2 + b)    *)

U2[-97., 177.]
(*    0.0000130837    *)

Such transformations exist very often for practical hypergeometric functions. Have a look at FunctionExpand and at Wolfram Research's massive database of hypergeometric transformations.