5
$\begingroup$

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 !

$\endgroup$
1
  • $\begingroup$ Or you can change the order of evaluation, e.g., Total@N[List @@ (f[10] // Apart)] evaluates to 3.02495 $\endgroup$
    – Bob Hanlon
    Apr 30, 2022 at 15:39

2 Answers 2

5
$\begingroup$

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}} 
$\endgroup$
5
$\begingroup$

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$.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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