# Getting wrong limit with Bessel

Bug introduced in 7.0 or earlier and 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])/Sqrt[1 + 24 n])/
(E^((Sqrt[n] π)/Sqrt)/(4 3^(1/4) n^(3/4))) -
(1 + (-((3 Sqrt)/(8 π)) + π/(48 Sqrt))/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)/(2048 π) - (35 π)/(36864 Sqrt) + π^3/(1990656 Sqrt) *)


This expression is equal to

(* 0.007709031447942101952246679 *)


But the correct result is

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


This is equal to

(* -0.009474863397773687074232327 *)


My question is: why is the term

(* -((315 Sqrt)/(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])/Sqrt[1 + 24 n])/
(E^((Sqrt[n] π)/Sqrt)/(4 3^(1/4) n^(3/4))) -
(1 + (-((3 Sqrt)/(8 π)) + π/(48 Sqrt))/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);


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])/
Sqrt[1 +
24 n])/(E^((Sqrt[n] π)/
Sqrt)/(4 3^(1/4) n^(3/
4))) - (1 + (-((3 Sqrt)/(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 π), (-1837080 + 102060 π^2 -
1890 π^4 + π^6)/(
1990656 Sqrt π^3), (-1837080 + 102060 π^2 -
1890 π^4 + π^6)/(
1990656 Sqrt π^3), (-1837080 + 102060 π^2 -
1890 π^4 + π^6)/(
1990656 Sqrt π^3), (-1837080 + 102060 π^2 -
1890 π^4 + π^6)/(1990656 Sqrt π^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}