A strange behaviour with function $\,_2F_1$

Question

I posted today a question on MSE 1

  u=(5*(49*Pi^2*Zeta[3] - 558*Zeta[5]))/(77*Pi^2*Zeta[3] - 930*Zeta[5])
  f[k_]:= Hypergeometric2F1[1/2, -k, 3/2, u]

Computing

  N[f[10]] =-8.29107
  N[f[10.]]=+3.02495

The second one is correct.

Could any one accept to make me understanging what is going wrong ?

Thanks in advance

Edit

After comments and answer, if I rationalize

 u=-137172947796693/512328299943181

the problem disappears.

Thanks, everybody !

Or you can change the order of evaluation, e.g., Total@N[List @@ (f[10] // Apart)] evaluates to 3.02495 — Bob Hanlon, Apr 30 '22 at 15:39

score 5 · Answer 1 · answered Apr 30 '22 at 14:47

It looks like numerical instability due to cancellations among terms in f[10] with the machine precision. If you specify the precision that you want to obtain, then it gives the correct answer.

N[f[10], 40]  (* 3.024945718826484765909050369309114702577 *)

Details: the sum in the numerator Numerator[f[10]] / 1835008 has 11 terms. You can see them by

List @@ (Numerator[f[10]] / 1835008)

or their numerical values with the machine precision:

N[List @@ (Numerator[f[10]] / 1835008)]

Summing up them gives a wrong result:

Plus @@ N[List @@ (Numerator[f[10]] / 1835008)]  (* -7.20576*^16 *)

From the following experiments, it seems that you need 19 digits or more to get the first digit correctly:

Table[{n,Plus @@ N[List @@ (Numerator[f[10]] / 1835008), n]},{n,16, 22}]

{{16,0``-17.31500855618509},{17,0.*^16},{18,3.*^16},{19,2.6*^16},{20,2.63*^16},{21,2.629*^16},{22,2.6290*^16}}

Michael E2 · Answer 2 · 2022-04-30T14:54:40.140

It is, sadly, simple round-off error. (I was hoping to learn something more interesting.)

Lets look at f[10]:

f10 = f[10]
(*    
  (1835008 (43551160253105097901 π^20 Zeta[3]^10 -
     4808402859534428003850 π^18 Zeta[3]^9 Zeta[5] +
     238905994054787308522350 π^16 Zeta[3]^8 Zeta[5]^2 -
     7034330305341815976769500 π^14 Zeta[3]^7 Zeta[5]^3 +
     135925793707604259244173750 π^12 Zeta[3]^6 Zeta[5]^4 -
     1801105572337397766462487500 π^10 Zeta[3]^5 Zeta[5]^5 +
     16574098514883746872868437500 π^8 Zeta[3]^4 Zeta[5]^6 -
     104587763224772152070746875000 π^6 Zeta[3]^3 Zeta[5]^7 +
     433131075279238064824511718750 π^4 Zeta[3]^2 Zeta[5]^8 -
     1063000622677773009307148437500 π^2 Zeta[3] Zeta[5]^9 +
     1174035206087269790337070312500 Zeta[5]^10)) /
  (138567 (77 π^2 Zeta[3] - 930 Zeta[5])^10)
*)

The terms being added/subtracted in the numerator and denominator are quite large:

{nf, df} = N@N[NumeratorDenominator[f10], $MachinePrecision]; 
terms1 = Numerator@f10 /. Plus -> List // N
terms2 = Denominator@f10 /. Plus -> List // N
(*
  {4.4145410^36, -4.2600110^37, 1.8499510^38, -4.760810^38, 
   8.0404910^38, -9.3120110^38, 7.4895910^38, -4.1307910^38, 
   1.4951810^38, -3.2072510^37, 3.0960210^36}
  {5.6080910^34, 9.6377210^34}
)

The expected round-off error is approximately:

1/df^2 (df*Abs@Max[terms1*$MachineEpsilon/2] Sqrt@Length@terms1 + 
   nf*Abs@Max[terms2*$MachineEpsilon/2] Sqrt@Length@terms2)
(*  18.5673  *)

Since f[10] is exact, one can use the adaptive-precision capability of N[] to compute an accurate value by giving it a numeric target precision in its second argument. For instance, we can use $MachinePrecision, which is 15.95... (but not the symbol MachinePrecision, which represents machine precision, the default target precision of N[expr]):

N@N[f10, $MachinePrecision]
(*  3.02495  *)

With f[10.] the internal code does not make the polynomial expansion, but goes straight to the numerical routine for calculating ${}_2F_1$.

A strange behaviour with function $\,_2F_1$

2 Answers2