No response to an infinite limit

Question

Trying to compute this limit, things seem to get frozen, even for hours. What would you
recommend to fastly compute it?

Limit[(Integrate[Product[x^k (1 - x^k), {k, 1, n}], {x, 0, 1}])^(1/n),n -> Infinity]

Answer 1

I very much suspect that the limit is not $1/4$, but rather $$\exp \left( \max_{y>0} \int_{t=0}^1 \left( -yt + \log(1-e^{-yt}) \right) dt \right) \approx 0.185155.$$

If I were writing this up on math.SE, I'd start right in on a proof of this, but on this site I think that it would be more welcome to show some of the Mathematica techniques I used to come up with this.

---

The first step in estimating an integral is to see where it is large.

integrand[x_, n_] := Product[x^k (1 - x^k), {k, 1, n}]
xMax[n_] := xMax[n] = x /. Last[NMaximize[{integrand[x, n], 0 < x < 1}, x]]
yMax[n_] := yMax[n] = First[NMaximize[{integrand[x, n], 0 < x < 1}, x]]
Do[{xMax[n], yMax[n]}, {n,5,30}]

Here xMax is where the function is largest, and yMax is its largest value. The NMaximize's take time, so I am computing the first bunch once and for all.

My first goal is to get an idea of the shape of the function. Right now, I don't care so much about its absolute size, I want to know what it looks like. Here is integrand, renormalized to take values from $0$ to $1$.

Plot[Table[integrand[x, 5*m]/yMax[5*m], {m, 1, 5}] , {x, 0, 1},PlotRange -> {0, 1}]

So we're talking about a spike near $1$, which gets narrower as $n \to \infty$. The width of that spike looks like $1/n^c$ for some small constant $c$ -- in any case not exponential -- so $\int_{x=0}^1 (\mbox{your function}) dx$ should be roughly $n^{-c} \max_{0 < x < 1} (\mbox{your function})$ and we should have $$\lim_{n \to \infty} \left( \int_{x=0}^1 (\mbox{your function}) dx \right)^{1/n} \approx \lim_{n \to \infty} \left( \max_{0 < x <1 } (\mbox{your function}) dx \right)^{1/n}.$$ This will be the first of many steps I won't make rigorous.

So, let's figure out how tall that spike is and, while we're at it, where it is. Here is a semi-log plot of the maximum value:

ListPlot[Table[Log[yMax[n]], {n, 5, 30}]]

Straight as a pin! So yMax is dropping exponentially, and the constant in that exponential is the constant you want. For the fun of it, let's do a best fit line:

Fit[Table[{n, Log[yMax[n]]}, {n, 5, 30}], {n, 1}, n]
(* Output 0.795824 - 1.6565 n *)

That $-1.66$ is pretty good: I'll later calculate that it is actually $-1.69$. For now, I'll continue the mathematical analysis.

The next step is to figure out where that maximum is. We saw above that xMax is approaching $1$, so the most useful number is 1-xMax[n]. Here is a log-log plot. (Disclaimer: I actually figured out what the answer should be with pencil and paper first, and did the plot to check my work.)

ListPlot[Table[{Log[n], Log[1 - xMax[n]]}, {n, 10, 30}]]

Looks linear again, suggesting that xMax is roughly $1-a n^{-b}$. Let's find out what $b$ is:

Fit[Table[{n, Log[1 - xMax[n]]}, {n, 10, 30}], {Log[n], 1}, n]
(* Output -0.0213984 - 0.908425 Log[n] *)

That $0.9$ is close enough to $1$ to suggest that, probably, xMax is acting like $1-a n^{-1}$. (I'm cheating: I already got an exponent of $1$ in my hand computation, so this is just confirming evidence.)

It is usually a good idea to make a change of variables which puts the maximum of the integrand in a position independent of $n$. Usually, I'd use $x=1-y/n$, but in this case $x$ is being raised to a lot of powers, so I went with $x=e^{-y/n}$. So our integrand is $$\prod_k e^{-k y/n} (1-e^{-ky/n}) = \exp \left( \sum_k - y (k/n) + \log(1-e^{-y k/n}) \right).$$ That sum looks like a Riemann sum, so (nonrigorously) $$\prod_k e^{-k y/n} (1-e^{-ky/n}) \approx \exp \left( n \int_{t=0}^1 \left(-yt + \log(1-e^{yt}) \right)dt \right).$$

How good is this approximation?

Table[
  Plot[
    {integrand[E^(-y/(5*m)), 5*m], 
     E^(5*m*NIntegrate[-y*t + Log[1 - E^(-y*t)], {t, 0, 1}])}, 
    {y, 0, 3}, PlotRange -> {0, yMax[5*m]}]  , {m, 1, 5}]

Not great. Here I show the $n=5$ and $n=25$ plots, omitting the intermediate ones. It looks like I missed a constant factor somewhere, or else the convergence is slow. But notice that it is only off by a factor of $2$ or $3$, even though the value of the function is changing by many orders of magnitude, and I nailed the basic shape of the function. Any constant factor will be swamped by the $1/n$ in the exponent when I take the limit.

In short, I believe that the integrand is approximately $\exp \left( n \int_{t=0}^1 \left(-yt + \log(1-e^{yt}) \right)dt \right)$. I also pick up a $(1/n) e^{-y/n} dy$ in place of $dx$, but that won't effect the asymtotics on this crude level. The largest value of the integrand should then be roughly $\exp \left( n \max_{y>0} \int_{t=0}^1 \left( -yt + \log(1-e^{-yt}) \right) dt \right)$ and the $n$-th root of the integral should approach $\exp \left( \max_{y>0} \int_{t=0}^1 \left( -yt + \log(1-e^{-yt}) \right) dt \right)$.

Before we compute that maximum, let's take a look at the function we are maximizing:

Plot[NIntegrate[-y*t + Log[1 - E^(-y*t)], {t, 0, 1}], {y, 0, 3}]

Good, nothing funny about it. A nice smooth peak. Incidently, this helps us analyze the width of the spike from back at the start: For $\exp \left( n \int_{t=0}^1 \left(-yt + \log(1-e^{yt}) \right)dt \right)$ to be above $1/e$ of it's maximum value, $\int_{t=0}^1 \left(-yt + \log(1-e^{yt}) \right) dt$ must be no less than $1/n$ below its maximum. Since that curve looks like it is smooth with a nonzero second derivative at it's max, this means that $y$ must be within $\pm 1/n^{1/2}$ of its optimum, so $x$ must be within $\pm 1/n^{3/2}$ of its optimum. In short, the width of the spikes from my first figure is probably $1/n^{3/2}$; a careful proof would have to establish this in order to justify approximating the integral by the maximum of its integrand.

There is only one thing left to do.

maxIntegral = NMaximize[{NIntegrate[-y*t + Log[1 - E^(-y*t)], {t, 0, 1}], 0 < y}, y]
(*  {-1.68656, {y -> 1.40505}} *)

E^First[maxIntegral]
(* 0.185155 *)

Answer 2

Even computing this numerically for large n is a bit of work..

integrand[x_?NumericQ, n_?IntegerQ] := Product[x^k (1 - x^k), {k, 1, n}]
ListPlot[Table[{n, 
       NIntegrate[integrand[x, n], {x, 0, 1}, WorkingPrecision -> 100, 
       MaxRecursion -> 20]^(1/n) }, {n, {1, 2, 3, 8, 10, 20, 50, 100, 1000}}], Joined -> True]

Looks to be approaching ~0.1841

Edit..perhpas useful.. We can series expand the product:

Series[ Product[x^k (1 - x^k), {k, 1, n}] , {n, Infinity, 1}]

This allows a somewhat faster numerical integration

n = 10000; 
NIntegrate[x^(n/2 + n^2/2) QPochhammer[x, x, n] , {x, 0, 1}, 
     WorkingPrecision -> 100, MaxRecursion -> 20]^(1/n)

-> 0.1850069

Also, since QPochhammer[x, x, large_n] is bounded between 0:1 we can do this:

 Limit[Simplify[ Integrate[x^(n/2 + n^2/2) (1)  , {x, 0, 1} ]^(1/n) , 
         Assumptions -> {n > 0}] , n -> Infinity] -> 1

Which I'd take as evidence that the original expression is convergent to something less than 1.

Answer 3

Just a note in the integration that takes a bit of CPU time. As the integrand is a polynomial its integral 0->1 can be obtained replacing $x^k$ by $1/(k+1)$ throughout. This will compute the integral by first computing the coefficient of the polynomial :

F[n_] := Module[{coefList},
   coefList = CoefficientList[Product[x^k (1 - x^k), {k, 1, n}], x];
   (Total[coefList[[#]]/# & /@ Range@Length@c[n] ])^(1/n) // N
   ];
F[20] // Timing 
(* {0.007211, 0.16846} *)
F[500] // Timing 
(* {78.636256, 0.183303} *)

whereas Integrate is way slower :

integrand[x_, n_] := Product[x^k (1 - x^k), {k, 1, n}]
With[{n = 20}, (Integrate[integrand[x, n], {x, 0, 1}])^(1/n) // N] // Timing
(* {9.616129, 0.16846} *)

F also appears to converge :

DiscretePlot[F[k], {k, 10, 200, 10}, Joined -> True]

BTW, the coefficients are rather neat, look at

   With[{n = 100},
     ListPlot[CoefficientList[Product[x^k (1 - x^k), {k, 1, n}], x], 
     PlotRange -> Full, ImageSize -> Large]]

Answer 4

Products are polynomials in x:

Table[Expand[Product[x^k*(1-x^k), {k, 1, n}]], {n, 1, 3}]
(* {x - x^2, x^3 - x^4 - x^5 + x^6, x^6 - x^7 - x^8 + x^10 + x^11 - x^12} *)

Integration of polynomials is easy:

Table[Integrate[Expand[Product[x^k*(1-x^k), {k, 1, n}]], x], {n, 1, 3}]
(* {x^2/2 - x^3/3, x^4/4 - x^5/5 - x^6/6 + x^7/7, x^7/7 - x^8/8 - x^9/9 + x^11/11 + x^12/12 - x^13/13} *)

For integral x=(0..1) set x=1:

n=1 1/2 - 1/3 = 1/6
n=2 1/4 - 1/5 - 1/6 + 1/7 = 11/420
n=3 1/7 - 1/8 - 1/9 + 1/11 + 1/12 - 1/13 = 293/72072
...

This efficient program generate a partial values of the integral (example is for the first six terms):

nmax=6; p=1; Table[p=Expand[p*x^n*(1-x^n)]; Total[CoefficientList[p,x]/Range[1,Exponent[p,x]+1]],{n,1,nmax}]
(* {1/6, 11/420, 293/72072, 487/760760, 129952159/1266697832400, 13084761625/783333734619744} *)

1000 terms took 23 minutes

nmax=1000; p=1; pint=Table[PrintTemporary[n]; p=Expand[p*x^n*(1-x^n)]; (Total[CoefficientList[p,x]/Range[1,Exponent[p,x]+1]]),{n,1,nmax}];
pintn=N[pint,200];

The convergence is slow

N[Table[pintn[[n]]^(1/n), {n, 995, 1000}], 40]
(* 
0.1840951809864800058349926417542938496702
0.1840960569507384522835520850517014404117
0.1840969313474495532267252012932767330364
0.1840978041809446369610699607582775721761
0.1840986754555388332987926055970487734842
0.1840995451755311499834383836140450183037
*)

Now I used the Richardson extrapolation and get a numerical values of the limit

funs[n_]:=pintn[[n]]^(1/n);
Do[Print[N[Sum[(-1)^(m+j)*funs[j*Floor[Length[pintn]/m]]*j^(m-1)/(j-1)!/(m-j)!,{j,1,m}],40]],{m,20,500,20}]
(*
0.1851478301101455130574734351538772362078
0.1851515730538468777493448434262869900998
0.1851525397534921153879509477106645265955
0.1851534344830868872547060203335950203739
0.1851538059763060513927789152414149494556
*)

Summary: I confirmed that $\lim\limits_{n\to\infty} \left(\int_0^1 \prod\limits_{k=1}^n x^k (1-x^k)\, \mathrm dx\right)^{1/n} = 0.18515\dots$

I added a new sequences to the OEIS, a values of the integral http://oeis.org/A258191 (numerators) http://oeis.org/A258192 (denominators)

---

The new results from 2000 terms (4 hours of CPU, 2 GB of RAM)

N[Table[pintn[[n]]^(1/n),{n,1995,2000}],40]
(*
0.1845615147056008261328531826856077053076
0.1845617654131831049853050447063095644066
0.1845620158931790428872902467157370250407
0.1845622661459067900133085427719476449285
0.1845625161716838953451706314205026095200
0.1845627659708273081040558246614154404272
*)

Richardson extrapolation:

funs[n_]:=pintn[[n]]^(1/n);
Do[Print[N[Sum[(-1)^(m+j)*funs[j*Floor[Length[pintn]/m]]*j^(m-1)/(j-1)!/(m-j)!,{j,1,m}],40]],{m,50,250,50}]
(*
0.1851515603160127958077453469352155726565
0.1851534313485228032278091602209464381079
0.1851540203809453934173356146638287154990
0.1851543619457164641835717384140678362809
0.1851545476774659555845959464451497858099
*)

Wynn epsilon:

SequenceLimit[pintn^(1/Range[Length[pintn]])]
(* 0.1851544656 *)

---

Conjecture: The limit is equal to

0.1851552893223595946473132111979542852738...
= 1 / 5.400871904118154152466091119104270052029...

see OEIS http://oeis.org/A107379

---

A similar limit: $\lim\limits_{n\to\infty} \left(\int_0^1 \prod\limits_{k=1}^n (1-x^k)\, \mathrm dx\right) = 0.368412535931433652321316597327851\dots$

after the Euler pentagonal theorem is the limit equal to

I created a new sequences http://oeis.org/A258229 (numerators), http://oeis.org/A258230 (denominators).

Answer 5

I don't know how to show the limit is 1/4. Ithink it can be bounded above by 1/4 as follows. Note that there is a gap in the argument below.

(1) Replace factors 1-x^k by 1-x^(n/2).

nn = Integrate[
  Product[x^k, {k, 1, n}]*Product[(1 - x^(n/2)), {k, 1, n}],
{x, 0, 1}, Assumptions -> n > 1000]

(* Out[147]= (2 Gamma[n] Gamma[1 + 2/n + n])/Gamma[2 (1 + 1/n + n)] *)

Limit[nn^(1/n), n -> Infinity]

(* Out[148]= 1/4 *)

(2) Investigate the error factor. (Also, check to make sure I didn't mess this up in some woeful manner.)

factor = 
 Limit[Product[(1 - x^(n/2))/(1 - x^k), {k, 1, n}], n -> Infinity, 
  Assumptions -> 0 < x < 1]

(* Out[152]= ((-1 + x)^2 (1 + x))/(QPochhammer[x, x] QPochhammer[x, x, 2]) *)

f2 = Together[FunctionExpand[factor, Assumptions -> 0 < x < 1]]

(* Out[176]= 1/QPochhammer[x, x] *)

(3) Show that this factor is bounded below by 1. This can be done by grouping the denominator in pairs of the form (1-x^k)*(1-x^(n-k)) and showing these are smaller than (1-x^(n/2))^2. Here I gotta run so this is a detail to check.

Answer 6

A better way how compute the maximum of integral

$$\int_{0}^1 \left( -yt + \log(1-e^{-yt}) \right) dt $$

is following. Derivative must be equal to zero at the maximum.

FullSimplify[D[Integrate[-y*t + Log[1 - Exp[-y*t]], t], y]]

$$\frac{\text{Li}_2\left(e^{t y}\right)+t y \left(\log \left(1-e^{t y}\right)-t y\right)}{y^2}$$

For t=1 we have

(t y (-t y + Log[1-E^(t y)]) + PolyLog[2,E^(t y)])/y^2/.t->1

$$\frac{\text{Li}_2\left(e^y\right)+y \left(\log \left(1-e^y\right)-y\right)}{y^2}$$

For t=0 we must use a limit

Limit[(t y (-t y + Log[1-E^(t y)]) + PolyLog[2,E^(t y)])/y^2, t->0]

$$\frac{\pi ^2}{6 y^2}$$

Maximum is at the point y, where y is the root of the equation $$-\text{Li}_2\left(e^y\right)+y^2-y \log \left(1-e^y\right)+\frac{\pi ^2}{6}=0$$

We have

FindRoot[Pi^2/6 + y^2 - y*Log[1-E^y] - PolyLog[2,E^y] == 0, {y,1}, WorkingPrecision->100]
(* {y->1.405054216578477533211873399082087670321440412190343393755288523553823298202008968431771181607706332} *)

The final value of the limit is

Re[Exp[Integrate[-y*t+Log[1-Exp[-y*t]] /.FindRoot[Pi^2/6 + y^2 - y*Log[1-E^y] - PolyLog[2,E^y]==0, {y,1}, WorkingPrecision->100], {t,0,1}]]]
(* 0.1851552893223595946473132111979542852738223644590750839856055303609697606866994281773438299262 *)

Or simplified

(r-1)/r^2 /.FindRoot[-PolyLog[2, 1-r] == Log[r]^2, {r, E}, WorkingPrecision->100]
(* 0.185155289322359594647313211197954285273822364459075083985605530360969760686699428177343829926166036 *)

$\lim\limits_{n\to\infty} \left(\int_0^1 \prod\limits_{k=1}^n x^k (1-x^k)\, \mathrm dx\right)^{1/n} = 0.18515528932235959464731321119795428527\dots$

This result agree with my previous conjecture, see http://oeis.org/A258234