# Why do these identical limits give different results?

Bug introduced in 12.0.0.0. and persisting through 13.1.0.0.

I wanted to calculate this limit:

 Limit[E^(E^ProductLog[n] (-1 + E^(1/2 Sqrt[ProductLog[n]/n]))) (1 + 1/2 Sqrt[1/(n*ProductLog[n])])^-n, n -> Infinity]


The correct result is E^(1/8)

Furthermore, it is true that

 Limit[(E^(-(1/2) E^(ProductLog[n]/2)))/ (1 + 1/2 Sqrt[1/(n*ProductLog[n])])^-n, n -> Infinity]


is equal to 1

Therefore, it must be

 Limit[E^(E^ProductLog[n] (-1 + E^(Sqrt[ProductLog[n]]/(2 Sqrt[n]))))*(E^(-(1/2) E^(ProductLog[n]/2))) , n -> Infinity]


also equal to E^(1/8) But the result shocked me. Mathematica gives Infinity !!!

The graph confirms my original assumption

 Plot[E^(E^ProductLog[n] (-1 + E^(Sqrt[ProductLog[n]]/(2 Sqrt[n]))))*(E^(-(1/2) E^(ProductLog[n]/2))) , {n, 1, 1000000}]


On deeper analysis, I was able to better isolate the problem and have the following simpler example.

 Limit[-(1/2) Sqrt[n/ProductLog[n]] + ((-1 + E^(1/(2 Sqrt[n/ProductLog[n]]))) n)/ ProductLog[n], n -> Infinity]


again gives the wrong result of Infinity

The limit is equal to 1/8, which is supported by the following graph

 Plot[-(1/2) Sqrt[n/ProductLog[n]] + ((-1 + E^(1/(2 Sqrt[n/ProductLog[n]]))) n)/ ProductLog[n], {n, 1, 1000000}]


• It depends on the version. The bug is in versions 12 and 13. Versions 10.2 and 11.0 calculate the limits correctly. Jul 8 at 20:04
• Please report to Wolfram Technical Support <support@wolfram.com>.
– Alan
Jul 8 at 20:25
• Do you have a question? It's not clear to me what you're asking for help with. Jul 8 at 21:06
• Question is in the title: Why do these identical limits give different results? Jul 8 at 21:18
• What do you mean, "why?" Do you want to know where the computation goes differently internally? That is, a comparison of the codes in the different versions? You should ask Wolfram about that. You seem to know the correct answer and you call it a bug in a comment. What more do you want to know? Jul 8 at 21:23

The easiest way, by a limited approach corresponding to approaching positive infinity:

expr = -(1/2) Sqrt[
n/ProductLog[n]] + ((-1 + E^(1/(2 Sqrt[n/ProductLog[n]]))) n)/
ProductLog[n];

Limit[expr /. n -> 1/n, n -> 0, Direction -> "FromAbove"]

(*  1/8  *)


Another workaround:

Limit[SeriesCoefficient[expr, {n, Infinity, 0}], n -> Infinity]

(*  1/8  *)


Possible cause of trouble:

The series at infinity seems divergent:

sd = Series[expr, {n, Infinity, 2}];
seriescoeffs = sd[[3]];
firstnumerator = sd[[4]];
denominator = sd[[6]];

kk = firstnumerator - 1;
(++kk; a[kk] (n^(kk/denominator)) ->
Limit[Simplify[# (n^(kk/denominator))],
n -> Infinity]) & /@ seriescoeffs
(*
{a[-1]/Sqrt[n] -> 0,
a[0]          -> 1/8,
Sqrt[n] a[1]  -> \[Infinity],
n a[2]        -> \[Infinity],
n^(3/2) a[3]  -> \[Infinity],
n^2 a[4]      -> \[Infinity]}
*)


A partial workaround.

\$Version

(* "13.1.0 for Mac OS X x86 (64-bit) (June 16, 2022)" *)

Clear["Global*"]

f[n_] = -(1/2) Sqrt[
n/ProductLog[n]] + ((-1 + E^(1/(2 Sqrt[n/ProductLog[n]]))) n)/
ProductLog[n] // FullSimplify

(* -(1/2) Sqrt[E^ProductLog[n]] +
E^ProductLog[n] (-1 + E^(1/(2 Sqrt[E^ProductLog[n]]))) *)


Using Asymptotic

asymp = Asymptotic[f[n], n -> Infinity] // Simplify[#, n > 1] &;

Limit[asymp, n -> Infinity]

(* 1/8 *)


Or using Series

ser = Series[f[n], {n, Infinity, 1}] //
Normal // Simplify[#, n > 1] &;

Limit[ser, n -> Infinity]

(* 1/8 *)

• Interesting that this gave me infinity: Series[expr, {n, Infinity, 1}] // Normal // Simplify // Limit[#, n -> Infinity] & b/c expr` wasn't simplified. Jul 8 at 23:25