12
$\begingroup$

Bug introduced in 7.0 or earlier and persisting through 11.1 or later. Fixed in 11.3


I computed a following limit (related to the asymptotic expansion of the sequence A000009 - number of partitions of n into distinct parts)

Expand[Limit[(((π BesselI[1, (Sqrt[1/24 + n] π)/Sqrt[3]])/Sqrt[1 + 24 n])/
             (E^((Sqrt[n] π)/Sqrt[3])/(4 3^(1/4) n^(3/4))) -
             (1 + (-((3 Sqrt[3])/(8 π)) + π/(48 Sqrt[3]))/Sqrt[n] +
              (-(5/128) - 45/(128 π^2) + π^2/13824)/n))*n^(3/2), n -> Infinity]]

Mathematica wrong output (in all versions from 7 to 11) is

(* (35 Sqrt[3])/(2048 π) - (35 π)/(36864 Sqrt[3]) + π^3/(1990656 Sqrt[3]) *)

This expression is equal to

(* 0.007709031447942101952246679 *)

But the correct result is

(* -((315 Sqrt[3])/(1024 π^3)) + (35 Sqrt[3])/(2048 π) -
   (35 π)/(36864 Sqrt[3]) + π^3/(1990656 Sqrt[3]) *)

This is equal to

(* -0.009474863397773687074232327 *)

My question is: why is the term

(* -((315 Sqrt[3])/(1024 π^3)) *)

missing in the Mathematica result? This is a bug!

Interesting is that numerically Mathematica evaluates a following expression correctly:

Table[N[(((π BesselI[1, (Sqrt[1/24 + n] π)/Sqrt[3]])/Sqrt[1 + 24 n])/
         (E^((Sqrt[n] π)/Sqrt[3])/(4 3^(1/4) n^(3/4))) -
         (1 + (-((3 Sqrt[3])/(8 π)) + π/(48 Sqrt[3]))/Sqrt[n] +
          (-(5/128) - 45/(128 π^2) + π^2/13824)/n))*n^(3/2), 10],
      {n, 1000000, 10000000, 1000000}]

(* {-0.009484688338, -0.009481807906, -0.009480532562, -0.009479772520, -0.009479253935, -0.009478871178, -0.009478573728, -0.009478333979, -0.009478135403, -0.009477967422} *)

Maple evaluates this limit correctly. Here is the Maple code:

expand(simplify(limit((Pi * BesselI(1, Pi * sqrt((n+1/24) * (1/3)))*(4 * 3^(1/4) * n^(3/4))/(sqrt(24 * n+1)*exp(Pi * sqrt((1/3) * n)))-1-(Pi/(48 * sqrt(3))-3 * sqrt(3)/(8 * Pi))/sqrt(n)-((1/13824) * Pi^2-5/128-45/(128 * Pi^2))/n) * n^(3/2), n = infinity))); evalf(%, 60);
$\endgroup$
  • $\begingroup$ Tested it anyone under new version 11.3 ? $\endgroup$ – Vaclav Kotesovec Mar 10 '18 at 22:17
12
$\begingroup$

I think it is something like a bug of Mathematica: It seems that Mathematica uses Series in order to calculate the Limit. When you take Series up to a expansion order that is high enough, you get the right answer. Mathematica doesn't take into account that the summation of Sqrt[1/n], n-> Infinity terms is not limited. When you take also the (1/n)3/2 term you get the right result. Look:

f[n_] = 
  (((π BesselI[1, (Sqrt[1/24 + n] π)/Sqrt[3]])/
    Sqrt[1 + 
      24 n])/(E^((Sqrt[n] π)/
       Sqrt[3])/(4 3^(1/4) n^(3/
         4))) - (1 + (-((3 Sqrt[3])/(8 π)) + π/(48 Sqrt[
          3]))/Sqrt[
     n] + (-(5/128) - 45/(128 π^2) + π^2/13824)/n))*
 n^(3/2) // Expand

Table[Series[f[n], {n, Infinity, i}], {i, 0, 3}]

It's to long to show here.

Table[Limit[Series[f[n], {n, Infinity, i}] // Normal, n -> Infinity], {i, 0, 6}]
{0, 0, (102060 - 1890 π^2 + π^4)/(
   1990656 Sqrt[3] π), (-1837080 + 102060 π^2 - 
   1890 π^4 + π^6)/(
  1990656 Sqrt[3] π^3), (-1837080 + 102060 π^2 - 
  1890 π^4 + π^6)/(
  1990656 Sqrt[3] π^3), (-1837080 + 102060 π^2 - 
  1890 π^4 + π^6)/(
   1990656 Sqrt[3] π^3), (-1837080 + 102060 π^2 - 
    1890 π^4 + π^6)/(1990656 Sqrt[3] π^3)}
Table[Limit[Series[f[n], {n, Infinity, i}] // Normal, n -> Infinity], {i, 0, 6}] // N
{0., 0., 0.00770903, -0.00947486, -0.00947486, -0.00947486, -0.00947486}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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