Limit of partial sums involving inverse squares

Question

Consider the finite sum

rs[x_, n_] := x/n Sum[n^2/(i + (n - i) x)^2, {i, 1, n}]

Is there a way to bring Mathematica to calculate the limit for n -> ∞?

I have tried Limit[] as well as NLimit[] without success.

Answer 1

For a fixed numerical value of x the sum limit can be found using Shanks transformation or Richardson extrapolation. The Richardson extrapolation transformation of the sequence gives faster convergence. The results seem to be in agreement with the plots in the answer by Willinski -- see the attached image.

Here is the code for Shanks:

Clear[Shanks]
Shanks[A_, n_] := (A[n + 2] A[n] - A[n + 1]^2)/(
  A[n + 2] + A[n] - 2 A[n + 1]);
Shanks[A_, 1, n_] := Shanks[A, n];
Shanks[A_, k_, n_] := (
  Shanks[A, -1 + k, n] Shanks[A, -1 + k, 2 + n] - 
   Shanks[A, -1 + k, 1 + n]^2)/(
  Shanks[A, -1 + k, n] + Shanks[A, -1 + k, 2 + n] - 
   2 Shanks[A, -1 + k, 1 + n]);

Here is the code for Richardson:

Clear[Richardson]
Richardson[A_, n_, N_] := 
 Total@Table[(A[n + k]*(n + k)^N *If[OddQ[k + N], -1, 1])/(
   k! (N - k)!), {k, 0, N}]

This code makes the table in the image:

rs[x_, n_] := x/n Sum[n^2/(i + (n - i) x)^2, {i, 1, n}]
sf[n_] := rs[1/2, n]

ns = {1, 2, 3, 4, 5, 6, 7, 15, 25, 50, 90};
res = Outer[#1[#2] &, {sf, Shanks[sf, #1] &, Shanks[sf, 3, #1] &, 
    Richardson[sf, #1, 1] &, Richardson[sf, #1, 6] &}, ns];
TableForm[Transpose[Map[N[#1, 20] &, res, {-1}]], 
 TableHeadings -> {ns, {"function value", "Shanks-1", "Shanks-3", 
    "Richardson-1", "Richardson-6"}}]

The book by Bender and Orszag "Advanced Mathematical Methods for Scientists and Engineers: Asymptotic Methods and Perturbation Theory" has good explanations of the Shanks and Richardson transformations. (See chapter 8.)

Answer 2

This post tackles the convergence acceleration of the Riemann integral in the same spirit as Anton's answer, except that I use a slight variation of one of the algorithms presented. In particular, I'm using this as an excuse to present the van den Broeck-Schwartz modification of the Wynn ε algorithm:

wgvs[seq_?VectorQ, h_: 1] := Module[{n = Length[seq], ep, v, w},
     Table[
           ep[k] = seq[[k]]; w = 0;
           Do[v = w; w = ep[j]; 
              ep[j] = v If[OddQ[k - j], h, 1] + 1/(ep[j + 1] - w),
              {j, k - 1, 1, -1}];
           ep[Mod[k, 2, 1]],
           {k, n}]]

The default setting of the second parameter corresponds to the classical Wynn algorithm.

For the OP's example:

rs[x_, n_] := x/n Sum[n^2/(i + (n - i) x)^2, {i, 1, n}, Method -> "Procedural"]
tab = Table[rs[1/2, 2^k], {k, 12}] // N;
res = wgvs[tab];

-Log10[Abs[tab - 1]]
   {0.51491, 0.770798, 1.04959, 1.33974, 1.6354, 1.93377, 2.23347, 2.53384,
    2.83454, 3.1354, 3.43635, 3.73734}

-Log10[Abs[res - 1]]
   {0.51491, 0.770798, 1.57696, 2.25408, 3.74516, 4.92739, 6.84944, 8.60221,
    10.9223, 13.2548, 15.6536, 14.9546}

where we see that Wynn ε achieved $\approx 14$ good digits with little additional effort.

For comparison, let's change the value of the second parameter of wgvs[] to 0; this corresponds to applying the iterated Aitken $\Delta^2$ process:

res = wgvs[tab, 0];
-Log10[Abs[res - 1]]
   {0.51491, 0.770798, 1.57696, 2.25408, 3.82643, 4.65844, 6.64954, 7.52831,
    8.61013, 10.7914, 13.1377, 14.8085}

---

In addition to Bender and Orszag, a good reference on convergence acceleration methods is Brezinski and Redivo-Zaglia's Extrapolation Methods: Theory and Practice. Weniger's survey paper is also a useful reference.

Answer 3

See EDIT #3 for a valid answer.

First answer

This is not an answer to my question but just a mathematical derivation of the limit the existence of which was even doubted in some comments and answers. I hope this does not spoil the creativity. The task for Mathematica is still open.

I came to the question considering something very elementary in calculus: the definition of the Riemann integral.

With this hint you surely can do the derivation by yourself. Here is mine:

The limit is equal to 1, as Willinski first calculated numerically.

Proof: consider the intergal

fi = Integrate[1/t^2, {t, 1, x}, Assumptions -> x > 1]

(* Out[395]= (-1 + x)/x *)

And now consider the definition of the Riemann integral of a function g[x] between x=a and x=b as a limit of this sum

gs[a_, b_, n_] := (b - a)/n Sum[g[a + (b - a) k/n], {k, 0, n - 1}]

Letting

{a -> 1, b -> x, g[t] -> 1/t^2}

we have

gs1[x_, n_] := (x - 1)/n Sum[1/(1 + (x - 1) k/n)^2, {k, 0, n - 1}]

This sum is equal to our sum except for the factor (x-1)/x (and the summand for k=n which, however, is negligible in the limit). Hence our limit is equal to x/(x-1) times fi, that is equal to 1. QED.

EDIT #1

One attempt to help Mathematica calculate the limit could be replacing the sum over i by an integral over i. Not strict but at least a heuristic approach.

Here we go

x/n Integrate[n^2/(i + (n - i) x)^2, {i, 1, n}, Assumptions -> {n > 1, x > 1}]

(*
Out[481]= ((-1 + n) x)/(1 + (-1 + n) x)
*)

Limit[%, n -> \[Infinity]]

(* Out[480]= 1 *)

We must then justifiy the replacement.

EDIT #2

Let us try to callculate symbolically the sum involved

s[n_, x_] := Sum[1/(i + (n - i) x)^2, {i, 0, n}]

Assuming[x > 1, Timing[s[n, x]]]

(* Out[27]= $Aborted *)

Mathematica won't do this for real x. Hence we adopt the strategy "guess from integers" which I use frequently in such cases : look for results for some interger values of x, and guess what the general expression might be.

Table[{x, Sum[1/(i + (n - i) x)^2, {i, 0, n}]}, {x, 1, 7}] // Column

$$\left( \begin{array}{cc} 2 & \psi ^{(1)}(-2 n)-\psi ^{(1)}(1-n) \\ 3 & \frac{1}{4} \left(\psi ^{(1)}\left(-\frac{1}{2} (3 n)\right)-\psi ^{(1)}\left(1-\frac{n}{2}\right)\right) \\ 4 & \frac{1}{9} \left(\psi ^{(1)}\left(-\frac{1}{3} (4 n)\right)-\psi ^{(1)}\left(1-\frac{n}{3}\right)\right) \\ 5 & \frac{1}{16} \left(\psi ^{(1)}\left(-\frac{1}{4} (5 n)\right)-\psi ^{(1)}\left(1-\frac{n}{4}\right)\right) \\ 6 & \frac{1}{25} \left(\psi ^{(1)}\left(-\frac{1}{5} (6 n)\right)-\psi ^{(1)}\left(1-\frac{n}{5}\right)\right) \\ 7 & \frac{1}{36} \left(\psi ^{(1)}\left(-\frac{1}{6} (7 n)\right)-\psi ^{(1)}\left(1-\frac{n}{6}\right)\right) \\ \end{array} \right)$$

Where $\psi ^{(1)}(z)$ = PolyGamma[1, z]

Our guess is

sg[n_, x_] := 1/(x - 1)^2 (-PolyGamma[1, 1 - n/(x - 1)] + PolyGamma[1, -((x n)/(x - 1))])

But still: no results for the limit

For real x > 1

Timing[Limit[n x sg[n, x], {n -> \[Infinity]}]]

(*
Out[26]= {211.21, {Limit[(
   n x (-PolyGamma[1, 1 - n/(-1 + x)] + 
      PolyGamma[1, -((n x)/(-1 + x))]))/(-1 + x)^2, n -> \[Infinity]]}}
*)

And even not for x = 2

Timing[Limit[n x sg[n, x] /. x -> 2, {n -> \[Infinity]}]]

(*
Out[28]= {18.1741, {Limit[2 n (-PolyGamma[1, 1 - n] + PolyGamma[1, -2 n]), 
   n -> \[Infinity]]}}
*)

Tough stuff for Mathematica!

EDIT #3

Now I've found a direct way to an answer my original question, viz. how we can bring Mathematica to calculate the limit.

This method generalizes to similar expressions of which we have to take the limit.

Letting

z -> x - 1;

our sum becomes

s := x/n Sum[1/(1 + k/n z)^2, {k, 0, n}]

We have to calculate the limit of s for n->\[Infinity].

Expanding the summand into a power series

1/(1 + k/n z)^2 -> Sum[Binomial[-2, m] k^m z^m/n^m, {m, 0, \[Infinity]}]

(*
Out[50]= 1/(1 + (k z)/n)^2 -> n^2/(n + k z)^2
*)

and exchanging the order of summation gives us

s1 := x/n Sum[Binomial[-2, m] z^m c[n, m], {m, 0, \[Infinity]}]

where

c[m_, n_] = Sum[k^m/n^m, {k, 0, n}]

(*
Out[53]= n^-m (0^m + HarmonicNumber[n, -m])
*)

Now we take the limit n->\[Infinity], taking into account the factor 1/n

clim = Limit[n^-1 c[m, n], n -> \[Infinity], Assumptions -> {m \[Element] Integers, m >= 0}]

(* Out[58]= 1/(1 + m) *)

and the m-sum becomes

s2 = Sum[x z^m Binomial[-2, m] clim, {m, 0, \[Infinity]}]

(* Out[59]= x/(1 + z) *)

Finally, replacing z gives

slim = s2 /. z -> x - 1

(* Out[61]= 1 *)

QED.

Answer 4

It's not the solution you are looking for and surely you have tried the same.

rs[x_, n_] := x/n Sum[n^2/(i + (n - i) x)^2, {i, 1, n}]
rs[1., \[Infinity]]

Infinity::indet: Indeterminate expression 0 [Infinity] encountered.

Sum::div: Sum does not converge.

Plot[rs[2., n], {n, 1, 100000}, PlotPoints -> 100]

Addition

<< NumericalCalculus`
NLimit[rs[2, n], n -> 3]

Infinity::indet: Indeterminate expression ComplexInfinity+ComplexInfinity encountered. >>

NLimit::notnum: The expression Indeterminate is not numerical at the point n == 4.`. >>

NLimit[2 n (PolyGamma[1, 1 - 2 n] - PolyGamma[1, 1 - n]), n -> 3]

Plot[PolyGamma[1, 1 - 2 n] - PolyGamma[1, 1 - n], {n, 1, 3}]

NLimit[PolyGamma[1, 1 - 2 n] - PolyGamma[1, 1 - n], n -> 3]

Infinity::indet: Indeterminate expression ComplexInfinity+ComplexInfinity encountered. >>

NLimit::notnum: The expression Indeterminate is not numerical at the point n == 4.`. >>

I would be surprised if there exists a limit.

We know, with Limit[] and NLimit[] we cannot obtain a solution. With DiscretePlot we can get an idea, but that is not a prove.

Table[DiscretePlot[rs[x, k], {k, 1000}, PlotRange -> All, 
  AxesOrigin -> {0, 0}], {x, {0.1, 1., 2., 5.}}]

I'm looking forward to the limiting procedure.

Answer 5

The answer

\begin{equation} \lim\limits_{n\rightarrow\infty}\frac{x}{n}\sum _{i=1}^n \frac{n^2}{(x (n-i)+i)^2}=\lim\limits_{a\rightarrow\infty}\left\{\begin{aligned} a^2 \left(\pi ^2 \csc ^2(\pi a x)-\log (a)-\pi \cot (\pi a)\right)&& x>1\\ 1&&x=1\\ -a^2 \log (a)&& 0\le x<1\\ a^2 \left(\pi ^2 \csc ^2(\pi a x)-\log (a)\right) && x<0 \end{aligned}\right. \end{equation} for $$a=\frac{n}{x-1}$$ if $$x\in\mathbb{Q}$$

Derivation

We first note that we can do the symbolic sum

Sum[1/(i + (n - i) 2)^2, {i, 1, n}]

PolyGamma[1, 1 - 2 n] - PolyGamma[1, 1 - n]

which is true for all $n\not\in\mathbb{Z}$. For integer $n$, we need to take the finite part; for example:

With[{n=5},
 Sum[1/(i+(n-i) 2)^2,{i,1,n}]
 ==Normal[Series[
    PolyGamma[1,1-2 k]-PolyGamma[1,1-k]
    ,{k,n,0}]]/.k->\[Infinity]
]

True

Since we are interested in $n\rightarrow\infty$ limit, we will not be bothered with $n\in\mathbb{Z}$ cases below, but they can be extracted as above.

As the next thing, we observe the following sequence:

Table[{x, Sum[1/(i + (n - i) x)^2, {i, 1, n}]}, {x, 1, 5}] // Column

{1,1/n}
{2,PolyGamma[1,1-2 n]-PolyGamma[1,1-n]}
{3,1/4 (PolyGamma[1,1-(3 n)/2]-PolyGamma[1,1-n/2])}
{4,1/9 (PolyGamma[1,1-(4 n)/3]-PolyGamma[1,1-n/3])}
{5,1/16 (PolyGamma[1,1-(5 n)/4]-PolyGamma[1,1-n/4])}

We can actually check this for non-integer $x\in\mathbb{Q}$ values too, they obey the observed pattern. I am not sure if we can enlarge the domain from $x\in\mathbb{Q}$ to $x\in\mathbb{R}$, hence I will simply assume $x\in\mathbb{Q}$ in the rest of the derivation: \begin{equation} \sum _{i=1}^n \frac{n^2}{(x (n-i)+i)^2}=\frac{n^2 \left(\psi ^{(1)}\left(1-\frac{x n}{x-1}\right)-\psi \left(1-\frac{n}{x-1}\right)\right)}{(x-1)^2} \end{equation}

Clearly, we can rewrite this equation in a more compact form for $a\equiv\frac{n}{x-1}$ which covers all cases for which the sum is divergent (x=1 case is the only convergent case at hand):

\begin{equation} \sum _{i=1}^n \frac{n^2}{(x (n-i)+i)^2}=a^2 (\psi ^{(1)}(1-a x)-\psi ^{(0)}(1-a)) \end{equation}

We can now consider four different cases: $x>1,1>x>0,x=0,0>x$. For each case, we series expand around $a\sim\pm\infty$ and take the leading piece. In the end, we get the answer written at the top of the post. As an example, consider one of the cases:

FullSimplify[Normal[Series[
  a^2 (-PolyGamma[0, 1 - a] + PolyGamma[1, 1 - a x])
  , {a, \[Infinity], -2}, Assumptions -> a > 0 && x > 1]],  a > 0 && x > 1]

$$a^2 \left(\pi ^2 \csc ^2(\pi a x)-\log (a)-\pi \cot (\pi a)\right)$$