# A strange behaviour with function $\,_2F_1$

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 ?

Edit

 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 Apr 30, 2022 at 15:39

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


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.41454*10^36, -4.26001*10^37, 1.84995*10^38, -4.7608*10^38, 8.04049*10^38, -9.31201*10^38, 7.48959*10^38, -4.13079*10^38, 1.49518*10^38, -3.20725*10^37, 3.09602*10^36} {5.60809*10^34, 9.63772*10^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$$.