Re: a tricky limit
- To: mathgroup at smc.vnet.net
- Subject: [mg15650] Re: [mg15327] a tricky limit
- From: Jurgen Tischer <jtischer at col2.telecom.com.co>
- Date: Sat, 30 Jan 1999 04:28:45 -0500 (EST)
- Organization: Universidad del Valle
- References: <199901080915.EAA03988@smc.vnet.net.>
- Sender: owner-wri-mathgroup at wolfram.com
Once more the tricky limit. I worked out a different solution to
calculate more digits. It's a bit involved so here it goes:
The problem was to solve Product[(1-x^p[n]/(n+1)),{n,2,inf}]/(1-x).
where p[n]=n-n/Divisors[n][[2]]. Taking Log it changes to the sum
-Log[1-x]+Sum[Log[1-x^p[n]/(n+1)],{n,2,inf}]. Now use the series for
Log[1-x] and Log[1-x^p[n]/(n+1)] to get
Sum[x^n/n,{n,inf}]-Sum[(x^p[k]/(k+1))^n/n,{n,inf},{k,2,inf}]. Now for
x<1 everything is absolutely convergent and wi can reorder to get:
x + x^2/2 -Sum[(x^p[k]/(k+1))^n / n,{k,2,inf},{n,2,inf} +
Sum[x^(n+1)/(n+1)-x^p[n]/(n+1),{n,2,inf}].
Let's start with one of the first sums, say for n==2. A lower bound at x
is given by
In[1]:= Sum[(x^k/(1 + k))^2/2, {k, 2, Infinity}]
Out[1]= (-4*x^2 - x^4 + 4*PolyLog[2, x^2])/(8*x^2)
In[2]:= Limit[%, x -> 1]
Out[2]= 1/24*(-15 + 2*Pi^2)
and an upper bound by
In[3]:= Sum[(x^(k/2)/(k + 1))^2/2, {k, 2, Infinity}]
Out[3]= (-4*x - x^2 + 4*PolyLog[2, x])/(8*x)
In[4]:= Limit[%, x -> 1]
Out[4]= 1/24*(-15 + 2*Pi^2)
So the limits are the same, and the same as Sum[1/k^2/2,{k,2,inf}], so
by monotone convergence we may interchange limits. Inspection of
higher n's (or knowledge of Riemann's Zeta function, see Help on Zeta),
leads to the formula
Limit[Sum[(x^p[k]/(k+1))^n/n ,{n,2,inf}],x->1] == Zeta[n]-1-2^-n
and we need to sum Sum[Zeta[n]-1 -2^-n,{n,2,inf}] Mathematica doesn't
want to solve it but in parts:
In[5]:= Sum[-(1/(2^n*n)), {n, 2, Infinity}]
Out[5]= 1/2*(1 - 2*Log[2])
In[6]:= Sum[(Zeta[n] - 1)/n, {n, 2, Infinity}]
Out[6]= 1 - EulerGamma
In[7]:= s1 = 1 + 1/2 - 1/2*(1 - 2*Log[2]) + 1-EulerGamma)//Simplify
Out[7]= EulerGamma+Log[2]
We are left with
Limit[Sum[x^(x+1)/(n+1),{n,2,inf}]-Sum[x^p[n]/(n+1),{n,48}],x->1].
The idea now is to extract from the sum the terms of one same second
divisor of n, similar to the sieve for finding primes. One problem is
that there is no closed formula for the sum over all terms with the
same divisor, one has to split the sum into sums of the type
Sum[x^(n(k-1)/k),{k,a,inf,b}], and even for small p this splits into
many sums. Now:
In[8]:= Sum[x^(n + 1)/(n + 1) - x^((n*(k - 1))/k)/(n + 1), {n, a,
Infinity, b}]
Out[8]= (x^(1 + a)*LerchPhi[x^b, 1, (1 + a)/b])/b - (x^((a*(-1 + k))/k)*
LerchPhi[(x^((-1 + k)/k))^b, 1, (1 + a)/b])/
b
Inspection of cases for different values of a and b shows that the
limit of the sum does not depend of a and is given by 1/b Log[1 - 1/k].
So the sum of all terms with second divisor k=Prime[m] is equal to
c/bLog[1-1/k], where b is the period for k, which is equal (without
proof, found by inspection) to the product \[CapitalPi] of all primes
up to k (included), and c is the number of multiples of k in [k,k+P-1]
with second divisor k. I calculated the first 8 c's and then searched
the sequence at Sloane's, (again without proof) c is equal to the
product of the primes-1 up to k-1. So we have for upper and lower
bounds ub, lb:
ub[n]: Sum[c[k]/b[k] Log[1-1/k],{k,n}]
lb[n]: ub[n]-Log[1-1/(n+1)]+Sum[c[k]/b[k]](Log[1-1/(n+1)])==
ub[n]-(1-Sum[c[k]/b[k]])Log[1-1/(n+1)]
where the meaning of the symbols is {n,p,s,s2}== {n, c[n]/b[n],
sum[c[k]/b[k],{k,n-1}], sum[c[k]/b[k] Log [1-1/k],{k,n-1}]}.
After summing up 106100002 terms (it took 45 hours) I found the
following bounds for the limit:
ub=2.29217369524, lb=2.29217369419, this is an error of 1.06 10^-9.
but I still don't know the 8th digit. At least my "best bet" 2.29217
coincides with this solution. Of course I was willing to spend another
(estimated) 45 hours to find this last digit, but there came a new
problem. Somewhere between 105 million and 106 million the speed went
down so calculation time of 100 000 terms rose from 100 seconds to 5000
seconds. This would mean some 100 days of calculating. (It seems that
the algorithm for Prime changes about that number.) So that's the end
for me for the moment.
Jurgen
If you happen to have a faster computer, you could run it on the
following:
{n,p,s,s2}={106100002,1.20410332625684268556672671364644320606271`30.243*^-11,
0.97388264577427072385607700884973105920738`30.3338,-0.\
440862266213089289557616316841531205845154`30.6699} li={};
While[True,{n,p,s,s2}=Nest[fu,{n,p,s,s2},100000];
AppendTo[li,{n,p,s,s2}];Print[{Exp[s1+s2],Exp[s1+s2+s
Log[1-1/(n+1)]]}]]
The Print is just so you know your calculation is still alive and if you
can stop.
Arnold Knopfmacher wrote:
>
> I wish to obtain a numerical estimate (say 8 decimal digits) of the
> limit as x tends to 1 from below of the function
> h[x]=(Product[(1-fm[x]/(m+1)),{m,2,Infinity}])/(1-x) where
> fm[x]=x^(m-m/d) and d is the smallest divisor of m that is greater than
> 1. The problem is that when I replace Infinity by say 1000 as the upper
> limit of the product, the function blows up near 1. Visual inspection
> of the graph of h[x] for 0<x<0.9 say, suggests that the limit should
> have a value around 2.1. Can anyone help?
>
> Thanks
> Arnold Knopfmacher
> Dept of Computational and Applied Math Witwatersrand University
> South Africa
- References:
- a tricky limit
- From: "Arnold Knopfmacher" <arnoldk@gauss.cam.wits.ac.za>
- a tricky limit