Possible bug in finite sum over inverse squares $\sum\limits_{i=1}^n \frac{1}{(x (n-i)+i)^2}$

Question

Revisiting the problem Limit of partial sums involving inverse squares I found another difficulty with Sum[]

Consider this sum

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

Here we assume x > 0, and n Integer > 0.

We have for example

s[2, 10]

(* Out[11]= 2920725891004177/54192375991353600 *)

But considering the symbolic evaluation gives

s[2, n]

(* Out[9]= PolyGamma[1, 1 - 2 n] - PolyGamma[1, 1 - n] *)

This is definitely a wrong result.

Numerically this becomes even more obvious:

Limit[%, n -> 10]

(* Out[10]= -∞ *)

I would consider this behaviour of Sum[] as a bug.

Answer 1

Note:

Please see the end of the post for more information. Since I still believe that it is not a bug, I did not change the main part of the post.

---

It is not a bug!

Unless you explicitly state otherwise, Sum evaluates generically, ignoring specific conditions. In your case, the result for generic $n$ is correct, and you can still get $n=10$ case by taking its finite part:

Normal[Series[
PolyGamma[1, 1 - 2 n] - PolyGamma[1, 1 - n]
, {n, 10, 0}]] /. n -> \[Infinity]

$\frac{2920725891004177}{54192375991353600}$

which is same as direct evaluation s[2, 10] with OP's command

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

If you prefer Sum not to compute generic cases (i.e. $n\in\mathbb{C}$) because you plan to use specific cases (i.e. $n\in\mathbb{Z}$), then you should set GenerateConditions -> True:

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

We now get OP's preferred behavior:

{sImproved[2, 10], sImproved[2, n]}

 {2920725891004177/54192375991353600, Sum[1/(i + 2 (-i + n))^2, {i, 1, n}, GenerateConditions -> True]}

EDIT:

I would like to clarify my points:

1. One should always set GenerateConditions to True to avoid results computed for generic situations. This is true for both Sum and Integrate!

2. Some integrations may not be doable directly if one sets GenerateConditions to True. If the problem is possible infinities, one may turn it off and try Regularization option which may help certain situations.For example

   Sum[x, {x, 1, \[Infinity]}, GenerateConditions -> True]
   

would not evaluate whereas

    Sum[x, {x, 1, \[Infinity]}, Regularization -> "Dirichlet"]

    (* -(1/12) *)

3. For analytic summations where neither GenereateCondition not Regularization is used, results may be incorrect for specific cases. My point is that this is not a bug!

Let us check the given sum in the question for $x=2$. It is

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

which Mathematica calculates if it is plainly given as

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

whereas Mathematica actually refuses to compute it if we specify above-mentioned options:

Sum[1/(2 (n - i) + i)^2, {i, 1, n}, GenerateConditions -> True]
Sum[1/(2 (n - i) + i)^2, {i, 1, n}, Regularization -> "Dirichlet"]

both remain unevaluated. I was earlier incorrect to claim that the result PolyGamma[1, 1 - 2 n] - PolyGamma[1, 1 - n] is generically correct! However, the fact that we can get the correct result from finite parts of this result is not a coincidence and I still take this result sufficiently fine, but mathematically speaking, it is not the correct analytical continuation. Nevertheless, we can understand why Mathematica choses this analytic continuation.

Let's first note that the correct analytic continuation (credit goes to OP) is given as

PolyGamma[1, n] - PolyGamma[1, 2 n] (*correct analytic continuation *)

which is actually related to Mathematica result

PolyGamma[1, 1 - 2 n] - PolyGamma[1, 1 - n] (*Naive Mathematica sum*)

as

FullSimplify[
 (PolyGamma[1, 1 - 2 n] - PolyGamma[1, 1 - n]) 
 == (PolyGamma[1, n] -   PolyGamma[1, 2 n]) 
     + \[Pi]^2 (-Csc[n \[Pi]]^2 + Csc[2 n \[Pi]]^2)
]
(*True*)

This explains why taking finite part of the result gives the expected solution and this is not specific to current situation at hand: I am only speaking from personal experience with Mathematica but summation/integration results containing additional singularities to be discarded is not something I first saw in this situation.

Of course, mathematically speaking, the result is simply incorrect at face value but we can see why Mathematica gives the result with all those additional poles. Consider following sequences of sums:

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

$$\frac{1}{(1-2 n)^2}+\frac{1}{4 (n-1)^2}\\ \frac{1}{(1-2 n)^2}+\frac{1}{(3-2 n)^2}+\frac{1}{4 (n-1)^2}\\ \frac{1}{4 (n-1)^2}+\frac{1}{(1-2 n)^2}+\frac{1}{(3-2 n)^2}+\frac{1}{4 (n-2)^2}$$

as the upper limit goes to infinity, we expect the resultant sum to have poles at all $\{n|2n\in\mathbb{Z}^+\}$. Hence it is only natural for mathematica to yield the analytic continuation

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

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

which has all those poles for half integer $n$. Again, this result is correct only generically as we should not have poles for $n>\frac{x-1}{x}m$ ($n>m/2$ in this case), but we did not ask Mathematica to generate conditions for which the result is valid, and we do need such poles for low enough $n$.

For the original sum Sum[1/(2 (n - i) + i)^2, {i, 1, n}], Mathematica simply sets $m=n$, which gives the correct analytic continuation plus those poles. As I showed in the original post, one can get the correct result by simply discarding those parts! And as we mentioned above, this follows from the equality:

FullSimplify[
 (PolyGamma[1, 1 - 2 n] - PolyGamma[1, 1 - n]) 
 == (PolyGamma[1, n] -   PolyGamma[1, 2 n]) 
     + \[Pi]^2 (-Csc[n \[Pi]]^2 + Csc[2 n \[Pi]]^2)
]
(*True*)

So given the situation, one may label the behavior as buggy. I beg the differ: Mathematica actually refuses to compute if we give it sufficient input:

Sum[1/(2 (n - i) + i)^2, {i, 1, n}, GenerateConditions -> True]
Sum[1/(2 (n - i) + i)^2, {i, 1, n}, Regularization -> "Dirichlet"]
Sum[1/(2 (n - i) + i)^2, {i, 1, n}, Assumptions -> n > 0]
(* none of them computes *)

Answer 2

Abstract

In addition to other answers I'd like to give a more analytically centered approach. I'll show that for $\operatorname{Re}(x')>0.$ it holds that $-(\psi^{(1)}(x'+n)-\psi^{(1)}(x'))=\sum_{i=1}^{n}\frac{1}{(i-1+x')^2}$. With that we can show why mathematica is "generically right" but wrong in the specific case considered here.

Proof

According to wikipedia it holds that

$$ -[\psi(x'+n)-\psi(x')]=-\sum_{i=0}^{n-1}\frac{1}{i+x'} $$

for $\operatorname{Re}(x')>0$. If we now take the derivative of left and right hand side, we obtain the polygamma functions found in the mathemtica expression on the left hand side.

$$ -[\psi^{(1)}(x'+n)-\psi^{(1)}(x')]=\sum_{i=0}^{n-1}\frac{1}{(i+x')^2} $$

For convenience we shift the expression on the right hand side a little and obtain the desired result.

$$ -[\psi^{(1)}(x'+n)-\psi^{(1)}(x')]=\sum_{i=1}^{n}\frac{1}{(i-1+x')^2} $$

Applying the result to the specific case

We set $x'=(1-2n)$ in the upper formula. This clearly violates the premise, because x' will be negative. However, we will proceed by evaluating the expressions formally to explain mathematicas result.

$$ \begin{split} -[\psi^{(1)}((1-n)-\psi^{(1)}(1-2n)]&=\psi^{(1)}(1-2n)-\psi^{(1)}((1-n)\\ &=\sum_{i=1}^{n}\frac{1}{(i-1+(1-2n))^2} \end{split} $$

The left hand side equals mathematica Polygammes[1,1-2n]-Polygamma[1,1-n] and the denominator of the sum equals (i + (n - i) x)^2 /. x -> 2. This can be checked with

(i - 1 + (1 - 2 n))^2 == (i + (n - i)*x)^2 /. x -> 2 // FullSimplify

which gives True.

To stress my point again $x'=-1+(1-2 n)$ is negative for any $n\geq0$. Since mathematica cannot know what x is going to be in advance I assume that it just uses some generic transformations, which do not apply in all cases and especially not in this case.

** Edit 1 **

I'm not quite sure if the strict requirement $\operatorname{Re}(x')>0$ is really necessary for the transformation to hold or if the weaker requirement $x'\neq0,-1,-2,-3,...$ is sufficient. The term $\psi^{(1)}(1-2n)-\psi^{(1)}((1-n)$ is nevertheless divergent at least for $n\in\mathbb{N}>0$. This means the result is indeed wrong.

** Edit 2 **

In this edit I will elaborate on the limiting procedure. I will suggest a limiting procedure different from the limiting procedure Limit[PolyGamma[1,1-2n] - PolyGamma[1,1-n], n -> 10] in the OP. The limit considered here, namely Limit[PolyGamma[1,1-z] - PolyGamma[1,1-z+n], z ->2n] with n=10, yields the same result as Sum[ 1/(i + (n - i) x)^2, {i, 1, n}] with n=10 and x=2. Furthermore, I will also explain why I deem the limiting procedure suggested in OP as unnatural.

Limit

I assume that mathematica analyzes the expression inside the sum, possibly by pattern matching, and then applies some generic transformation. The sum can also be written as

$$ \tag{1}\label{sum-equation} \sum_{j=1}^{n}\frac{1}{(i+(n-i)x)^2}=\sum_{j=1}^{n}\frac{1}{(ia+z)^2} $$ with $a=1-x$ and $z=nx$.By evaluating the right hand side Sum[1/(j*a + z)^2, {j,1,n}], we obtain

$$\tag{2}\label{sum-polygamma-relation} \sum_{j=1}^{n}\frac{1}{ia+z}=\frac{1}{a^2}\left(\psi^{(1)}(1+z/a)-\psi^{(1)}(1+z/a+n)\right) $$ I will now take the right hand side and treat $z$ as the idependent variable rather than $n$. The left hand side has poles at $z=-\frac{1}{a},-\frac{2}{a},...,-\frac{n}{a}$. The right hand side has those poles as well, but additionally it has points where it is not defined. Those points are at $z=-\frac{j}{n}$ with $j\in\mathbb{N}>n$. Now if we want to visualize that we pick x=2 and n=10 this gives us a=-1 and z=-20. We can now plot either the righthandside or lefthandside

lefthandside = Sum[1/(-i + z)^2, {i, 1, n}];
righthandside = PolyGamma[1, 1 - z] - PolyGamma[1, 1 - z + n];
Show[Plot[lefthandside, {z, -22, 15}, PlotRange -> {0, 15}, AspectRatio -> 1], Graphics[{Red, Circle[{-20, 0.0539}, 0.5]}], PlotRangePadding -> {{0, 0}, {1.5, 0}}, AspectRatio -> Automatic]

and arrive at an identical plot for both. The red circle marks left- and right hand sides value at z=-20 which is about 0.0539. Note that the right hand side of \eqref{sum-polygamma-relation} is not defined at $z=-1,-2,-3,...$. Nevertheless it becomes obvious from the plot that there exists a continuous extension to all of $\mathbb{C}$. This continuous extension is given by the sum. When excluding the singular points left and right hand side agree on all of $\mathbb{C}$.

Conclusion

IMHO considering $n$ as the independent variable in the limiting procedure seems unnatural as non-integer, negativ or complex values for $n$ do not make sense as the upper limit of the sum. Hence it seems more natural to rewrite the sum in terms of polygamma functions as suggested in \eqref{sum-polygamma-relation} and evaluate those in a neighborhood of $z=xn$. For $x\not\in\mathbb{N}$ the polygamma functions given by mathematica agree with the sum. For $x\in\mathbb{N}$ the right hand side of \eqref{sum-polygamma-relation} is ill defined. Contrary to my previous assertion in "edit 1" the polygamma term does not diverge there. It is just not defined. Since the expression given by mathematica is only wrong on a subset of $\mathbb{C}$ one may or may not deem the mathemtica solution as "generically" correct.

Answer 3

Abstract

In this answer I'm going to use a numerical approach to answer the question and convince graphically that the sum and the difference of the polygamma functions agree in all of $\mathbb{C}$.

Intro

As stated in my other answer I assume that $-(\psi^{(1)}(z+n)-\psi^{(1)}(z))=\sum_{i=1}^{n}\frac{1}{(i-1+z)^2}$ holds for some $z\in\mathbb{C}$.

To get a rough idea how the left and right hand side look we plot both in the complex plane

ComplexPlot[-(PolyGamma[1, z + n] - PolyGamma[1, z]), {z, -30 - 10*I, 10 + 10*I}]

ComplexPlot[Sum[1/(i - 1 + z)^2, {i, 1, n}], {z, -30 - 10*I, 10 + 10*I}]

To my suprise, they look like they agree in the Subset of $\mathbb{C}$ we are interested in.

What about the divergence?

The sum gives

n=10;
Sum[1/(i - 1 + z)^2, {i, 1, n}] /. z -> (1 - 2*n) // N

0.0538955

The difference of the two polygamma functions gives near the point we are interested in:

c = 10^-8;
-(PolyGamma[1, z + n] - PolyGamma[1, z]) /. z -> (1 - 2*n) + c //N

0.0538955

I'm fairly confident now that in the limit $c\to0$ the divergent parts in the polygamma functions cancel.

Answer 4

This self answer consists of three parts. In the first past we provide closed expression for the sum for all vaues of the parameter $x \gt 0$, the second part shows a close relation between the correct result and the wrong one returned by Mathematica, finally, I briefly discuss a well appreciated but wrong answer.

EDIT 19.05.20: official rejection of bug statement from Wolfram added.

Closed expressions for the sum

To begin with, let me write down these closed expressions for the sum

si[x_?# > 1 &, n_] = (
      PolyGamma[1, n/(x - 1)] - PolyGamma[1, (n x)/(x - 1)])/(-1 + x)^2; (1)

si[x_?# < 1 &, n_] = (-PolyGamma[1, 1 + n/(1 - x)] + 
       PolyGamma[1, 1 + (n x)/(1 - x)])/(-1 + x)^2; (2)

and, of course,

si[1, n_] = 1/n^2; (3)

Checking numerically,

Table[{s[1/2, n], si[1/2, n]}, {n, 1, 5}]

(* Out[223]= {{1, 1}, {25/36, 25/36}, {469/900, 469/900}, {73249/176400, 
  73249/176400}, {547129/1587600, 547129/1587600}} *)

Table[{s[2, n], si[2, n]}, {n, 1, 5}]

(* Out[224]= {{1, 1}, {13/36, 13/36}, {769/3600, 769/3600}, {26581/
  176400, 26581/176400}, {737641/6350400, 737641/6350400}} *)

Now for x=2 we get

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

This is completely well defined and gives correct numerical results for n>0, in sharp contrast to the wrong generic result of Mathematica provided in the OP

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

Derivation

We start with deriving an integral representation of the sum.

Using the formula

Integrate[Log[1/z] z^(q - 1), {z, 0, 1}, Assumptions -> q > 0]

(* Out[77]= 1/q^2 *)

to replace the denominator, and doing the sum under the integral leads to

si0[x_, n_] := Integrate[Log[1/z] (z^n - z^(n x))/(z - z^x), {z, 0, 1}]

Now adding the assumption x>0 to the integral gives

Integrate[Log[1/z] (z^n - z^(n x))/(z - z^x), {z, 0, 1}, 
 Assumptions -> x > 0]

(* Out[146]= ConditionalExpression[(
 PolyGamma[1, n/(-1 + x)] - PolyGamma[1, (n x)/(-1 + x)])/(-1 + x)^2, 
 x > 1] *)

which returns the condition x>1 and the (correct) expression si[x>1,n] above.

But what happens in the remaining region of $x$?

Integrate[Log[1/z] (z^n - z^(n x))/(z - z^x), {z, 0, 1}, 
 Assumptions -> 0 < x < 1]

(* Out (-PolyGamma[1, n/(-1 + x)] + PolyGamma[1, (n x)/(-1 + x)])/(-1 + x)^2 *)

At $x=1/2$ this gives

4 (-PolyGamma[1, -2 n] + PolyGamma[1, -n])

which is just as wrong as the generic result sg[x,n].

Hence also Integrate[] has a similar bug as Sum[].

What's the remedy? Well, one idea is to just set x->1/2 before doing the integral:

Integrate[
 Log[1/z] (z^n - z^(n x))/(z - z^x) /. x -> 1/2, {z, 0, 1}, 
 Assumptions -> x > 0]

(* Out[219]= ConditionalExpression[
 4 (PolyGamma[1, 1 + n] - PolyGamma[1, 1 + 2 n]), Re[n] >= -(1/2)] *)

Now the result corresponds to the general correct formula si[x,n].

In order to find a general expression for the region 0<x<1 we could try several values of x and then guess the result. But here's a simpler way: let x->pi /4, and after evaluation replace back pi -> 4 x.

This heuristic trick gives si[x,n] valid for 0<x<1.

In summary: there is a bug (at least in my version 10.1.0) which causes wrong results the generic Sum[] as well as in the generic Integrate[].

But nevertheless we were able to find exact results using Mathematica, giving it a nudge.

Relation between correct and wrong expression

Let us compare these two expressions for $x = 2$

s0 = PolyGamma[1, 1 - 2 n] - PolyGamma[1, 1 - n] (* wrong *);
s1 = PolyGamma[1, n] - PolyGamma[1, 2 n] (* correct *);

s0 is returned by Mathematica, and it is wrong, as it leads to infinite values for positive integer n, s1 was derived by analytic continuation, and it is true as it gives the correct values for positive integers n.

For brevity we shall say s0 is wrong if it does not conincide with s1.

It is easy to show that s0 is wrong not only for positive integers n but also for real values, say n=1/5.

{s0, s1} /. n -> 1/5
% // N

(* Out[350]= {PolyGamma[1, 3/5] - PolyGamma[1, 4/5], PolyGamma[1, 1/5] - PolyGamma[1, 2/5]}

Out[351]= {1.33674, 18.992}
*)

So the problem can be decoupled from the divergence.

But we can derive the difference quite generally using the reflection formula of the polygamma function (https://de.wikipedia.org/wiki/Polygammafunktion)

$$(-1)^m \Psi_m(1-z) = \Psi_m(z) -\pi \frac{d^m}{dz^m}{\cot(\pi z)$$

which gives

s1 - s0 = \[Pi]^2/Sin[n \[Pi]]^2 - \[Pi]^2/Sin[2 n \[Pi]]^2

This quantity vanishes only for n=k[PlusMinus]1/3, k[Element]Integers.

Hence s0 is incorrect for all real n except countably many values. This is sometimes stated as "almost always inccorect".

Comment on the strange answer of Soner

I make this comment because this answer was pretty much appreciated by the community.

In spite of this support this answer is planily wrong, and this goes for the result as well as for the method.

The answer starts with the bold face line "it is not a bug", and then states: "Unless you explicitly state otherwise, Sum evaluates generically, ignoring specific conditions. In your case, the result for generic n is correct" and he quotes s0.

In fact the expression s0 is wrong for any positive integer n. Simply because it is divergent. And we have furthermore shown that s0 is wrong almost everywhere.

But let's nevertheless turn to the "reasoning" of Soner:

I'm not a newbee but I must admit that I haven't seen yet a magic procedure like this in oder to calculate the value of a function.

The first step is clear: an expansion in n about the value n=10 in question. To first order this yields

Normal[
 Series[PolyGamma[1, 1 - 2 n] - PolyGamma[1, 1 - n], {n, 10, 0}]]

(* Out[251]= 2920725891004177/54192375991353600 - 3/(4 (-10 + n)^2) *)

This shows clearly that the function has a double pole at n=10 and hence is divergent in this point.

Now the magic (or was it perhaps meant as a test of the audience?): the author just throws away the divergent part (by suddenly applying an additional limit n->inf, remember that we were at n=10) and declare the rest as the value of the function.

Simply stated: this is not mathematics. With the same "method" he could as well "prove" that Zeta[1] = 0.

And, of course, this answer cannot serve to reject the bug statment.

Official statement of Wolfram

Here is the official answer to my bug report

Betreff: Re: [CASE:4544347] Bug report Sum[], Integrate[]
Datum: Mon, 18 May 2020 12:09:25 -0500
Von: Wolfram Technical Support

"Hello Wolfgang,

This behavior is not a bug. As is documented in several places in the documentation

https://support.wolfram.com/39071?src=mathematica https://reference.wolfram.com/language/ref/FullSimplify.html#482986235 and https://reference.wolfram.com/language/ref/Sum.html#87823560

Mathematica assumes that unspecified symbolic variables are in general complex. Your symbolic expression fails at integers, but this is a measure-zero set of the complex plane, so the 'generic' validity of the Sum is upheld.

Further information on how to use GenerateConditions and Assumptions to avoid these issues is also available in the documentation at the included links."

My comment: To me this answer is remarkable but neither satisfactory nor useful for practical purposes. Considering that I have provided a Mathematica expression which yields the correct values also at integers the simple question is: why doesn't Mathematica return this expression? (Also, the Integrate question was not answered.)

Bottom line: I have received two different strange answers to my observation.

Answer 5

Modify the input under consideration:

Limit[s[2, n], n -> 10.]

Making the value of n in the limit Reals give a slightly different and much more instructive output.

Mathematica assigns internally a jump at the Integers if the sum invaried in the second argument:

Plot[s[2, n], {n, 9.5, 10.5}]

So for n smaller than 10 the sum and s[x,9] is in output and for n bigger than 10 the sum s[x,10] is in output. Both have different values

{s[2, 9], s[2, 10]}

(* {9064791593257/150117385017600, 2920725891004177/54192375991353600} *)

This is a finite jump as the Graphics shows, but in the output of the Limit built-in there is infinity as a replacement for indefinite, undefined. If used over Reals that are more appropriate for the handling of the second argument in use of Limit, than the output is more reliable as shown in the screenshot.

Mathematica uses a standard definition set for such functions as PolyGamma stemming from NIST most probable. It is always a best practice workaround to look up the definition on the corresponding NIST definition pages:

[PolyGamma

As can be read off the screen of that page, PolyGamma is not defined on Integers, because the denominator series shows a singularity for finite series members.

If that is not a satisfaction have a look at PolyGamma at NIST.

The convention depends strong on the use of the definition. There are references in chap. 5.21 for the computation of the Gamma and PolyGamma function including surveys. Another look worth is the book: NIST Handbook of Mathematical Functions Hardback and CD-ROM. Suitable might be Chap. 5.21 and 5.22. Or the book Nancy Blachman, Mathematica Quick Reference, Version 2

This is not just a problem around with Mathematica other numerical programs do alike.