Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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]
share|improve this question
add comment

4 Answers

up vote 8 down vote accepted

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}]

enter image description here

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}]]

enter image description here

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}]]

enter image description here

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}]

enter image description here enter image description here

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}]

enter image description here

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 *)
share|improve this answer
    
So (BesselI[0, 1] + Log@2)^Catalan? –  belisarius Dec 17 '13 at 4:05
    
Very nice. And an upvote (I mention this in part to surpass the 15 character minimum). –  Daniel Lichtblau Dec 17 '13 at 15:26
1  
@belisarius How did you come to that? –  ziyuang Dec 17 '13 at 20:52
    
@ziyuang http://oldweb.cecm.sfu.ca/cgi-bin/isc/lookup?number=0.185155&lookup_type=simple‌​ –  belisarius Dec 24 '13 at 0:17
add comment

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.

share|improve this answer
    
Yeah, it seems the limit exceeds $1/(2e)$. +1 for showing that. Finally, the limit should hit $1/4$ as expected. –  Chris's sis Dec 16 '13 at 20:21
add comment

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.

share|improve this answer
1  
We can use AM-GM inequality to show the upper bound $$x^k (1-x^k) \le \left(\frac{x^k+1-x^k}{2}\right)^2=\frac{1}{2^2}$$ –  Chris's sis Dec 17 '13 at 8:39
add comment

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]]

(I updated this plot with <code>PlotRange->Full</code>

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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