10
$\begingroup$

EDIT 11.10.2017

In order to avoid confusion which might arise from the statements I made in this question I declare here that the extensive discussion here has shown that the described behaviour is not a bug in Mathematica but it is a surprising mathematical fact which we have come here to call "miracle".

It might nevertheless be interesting to follow the pretty tenacious discussion in which the community guided me eventually to this conclusion.

Finally, I'd like to point out that to understand and explain this "miracle" and the accompanied "magic" number 7 is no longer an open question (see my self-answer 4).

Original post

I find that Mathematica returns a wrong result for sums of the type

$$s(n)=\sum_{k=0}^\infty \left(\frac{\sin (k)}{k}\right)^n$$

for n = 7, 8, ..., and suspect a bug here.

$Version

(* Out[38]= "10.1.0 for Microsoft Windows (64-bit) (March 24, 2015)" *)

Let

f[k_]= Sin[k]/k 

and consider the infinite sum

s[n_] := Sum[f[k]^n, {k, 0, ∞}]

Mathematica calculates correctly the first 6 values

Table[{n, s[n]}, {n, 1, 6}] // Expand

(* Out[15]= {{1, 1/2 + π/2}, {2, 1/2 + π/2}, {3, 
  1/2 + (3 π)/8}, {4, 1/2 + π/3}, {5, 
  1/2 + (115 π)/384}, {6, 1/2 + (11 π)/40}} *)

but it fails at n = 7, and it takes appreciably longer time to come up with

s[7]
N[%, 10]

(* Out[22]= 1 + (-23040 + 129423 π - 201684 π^2 + 
  144060 π^3 - 54880 π^4 + 11760 π^5 - 1344 π^6 + 
  64 π^7)/46080 *)

(* Out[23]= 1.302724351 *)

Which obviously has a completely different symbolic structure.

The correct result can be computed in Mathematica via the corresponding integral and making use of the Euer-Maclaurin formula which tells us in this case that the difference beween sum and integral is exactly 1/2 to give

1/2 + Integrate[f[k]^7, {k, 0, ∞}]
N[%, 10]

(* Out[25]= 1/2 + (5887 π)/23040 *)

(* Out[26]= 1.302715102 *)

An error appears also for n=8..10 so that I conclude that Mathamtica cannot deal correctly with the sum from n=7 onwards.

For the These values of n the wrong sum is greater that the correct result.

Added remark

I believe that there is some switch in the calculation method of Sum[] starting at at n=7, indicated also by the longer runtime mentioned above.

This is confirmed in Wolfram Alpha where the cases n = 1..6 are returned in almost no time whereas n=7 needs excess time which in turn requires Alpha Pro.

$\endgroup$
  • 1
    $\begingroup$ MMA 11.2 give the same symbolic answer what MMA 10.1.0.It seems not a bug.Check :1 + NSum[f[k]^7, {k, 1, \[Infinity]}, Method -> "AlternatingSigns", NSumTerms -> 10^5, WorkingPrecision -> 40].Maybe Euer-Maclaurin formula not working for n=7 ? $\endgroup$ – Mariusz Iwaniuk Oct 9 '17 at 20:03
  • $\begingroup$ @Mariusz Iwaniuk I gather from your remark, that the bug also appears in MMA 11.2. And IMHO, I would prefer to accept a bug in MMA than in the Euler-Maclaurin formula ;-) $\endgroup$ – Dr. Wolfgang Hintze Oct 9 '17 at 20:27
  • $\begingroup$ Here's what I get.. It looks to me like the sequence of partial sums converges to the value that the infinite sum spits out. $\endgroup$ – march Oct 9 '17 at 21:17
  • $\begingroup$ @march Isn't that a self-referential calculation? $\endgroup$ – Dr. Wolfgang Hintze Oct 9 '17 at 21:30
  • $\begingroup$ I don't think so. The methods used behind the scenes are different for infinite sums versus finite sums, since for finite sums Mathematica can just add the values, whereas there has to be some analytic machinery behind the scenes in order to compute the infinite sum. $\endgroup$ – march Oct 9 '17 at 21:39
7
$\begingroup$

Fourth answer

Introduction

Carl Woll pointed out in his solution that the series expansion around x = 0 is wrong for n=7. Alas, I have to confirm this and show it below. Hence my "proof" given in the answers 2 and 3 is incorrect, and probably its conclusion as well.

Hence it seems that the tide has swept back to my first answer.

As a little compensation for the to and fro the constructive part of this 4th answer will be a closer look at the cause of the "miracle".

This closer look paid, and I'm happy to present the explanation of the "magic seven" below.

Calculations

We restart with the x-dependent sum

s7a = Sum[(Sin[k x]/k)^7, {k, 1, \[Infinity]}]

(* Out[148]= 1/128 I (35 PolyLog[7, E^(-I x)] - 35 PolyLog[7, E^(I x)] - 
   21 PolyLog[7, E^(-3 I x)] + 21 PolyLog[7, E^(3 I x)] + 
   7 PolyLog[7, E^(-5 I x)] - 7 PolyLog[7, E^(5 I x)] - 
   PolyLog[7, E^(-7 I x)] + PolyLog[7, E^(7 I x)]) *)

Its value at x=1 is

N[s7a /. x -> 1, 20]

(* Out[160]= 0.30272435072915288381 + 0.*10^-21 I *)

Whereas my value was different

N[5887 \[Pi] /23040 - 1/2, 20]

(* Out[161]= 0.30271510206957954839 *)

This difference is due to the non analyticity. This can be seen in the 6th derivative

s7ap = D[s7a, {x, 6}]

(* Out[152]= 1/128 I (35 Log[1 - E^(-I x)] - 35 Log[1 - E^(I x)] - 
   15309 Log[1 - E^(-3 I x)] + 15309 Log[1 - E^(3 I x)] + 
   109375 Log[1 - E^(-5 I x)] - 109375 Log[1 - E^(5 I x)] - 
   117649 Log[1 - E^(-7 I x)] + 117649 Log[1 - E^(7 I x)]) *)

enter image description here

We can see that the point x=1 is above the jump of the derivative which is located at x~=0.897 and has a size of about 2402. Hence the point x=1 cannot be reached from a series expansion about x=0.

Comparing the cases n=6 and n=7 in one plot shows the essential difference: for n=6 (and lower values not shown here) the jump is above x=1 whereas for n=7 (and higher values) the jump happens below x=1.

In other words: for n<=6 the point x=1 can be reached from a series expansion around x=0, for n>=7 not.

enter image description here

As a rather extreme case consider n=15

enter image description here

From this we can understand at least part of the "miracle" of the magic number 7.

Is is worthwhile to mention that the "miracle" seems to be a "quantum effect", as it does not apear in the continuous integral version.

Explanation of the magic seven

Continuing the study of the jumps in the (n-1)st derivative I found an explanation for the magic seven which can be summarized best with the comment of bbgodfrey who stated yesterday dryly:

"I believe that 7 is magic, because it is the first integer greater than 2 Pi."

(this was more than a belief, I believe).

In other words

"7 is magic because 2[Pi]/7<1<2[Pi]/6"

I have already seen in the graphs that at n=7 the jump occurs for the first time below x = 1.

Let us define this situation as the "magic moment" or the "miracle".

We shall now make this quantitave.

Let us take a closer look at the (n-1)st derrivative of the symmetric sums for n=6 and n=7.

First n=6

ss6xp = D[Sum[(Sin[k x]/k)^6, {k, -\[Infinity], \[Infinity]}], {x, 5}]

(* Out[760]= (1/6048)(4354560 x - 90720 I Log[1 - E^(-2 I x)] + 
  90720 I Log[1 - E^(2 I x)] + 1161216 I Log[1 - E^(-4 I x)] - 
  1161216 I Log[1 - E^(4 I x)] - 1469664 I Log[1 - E^(-6 I x)] + 
  1469664 I Log[1 - E^(6 I x)]) *)

This function has branch points in x where the argument of the log-function is zero. This happens the first time for x>0 when the exponent is an integer multiple of 2[Pi]. Since we have 3 summands in pairs we find three branch points $\beta$ at

\[Beta]6 = 2 \[Pi]/{2, 4, 6}
% // N

(* Out[690]= {\[Pi], \[Pi]/2, \[Pi]/3} *)

(* Out[691]= {3.14159, 1.5708, 1.0472} *)

Notice that all branch point positions are greater than 1.

Now for n=7 we have in the same manner

ss7xp = D[Sum[(Sin[k x]/k)^7, {k, -\[Infinity], \[Infinity]}], {x, 6}]

(* Out[761]= -5040 x + 2 (5040 x + 
    1/128 I (35 Log[1 - E^(-I x)] - 35 Log[1 - E^(I x)] - 
       15309 Log[1 - E^(-3 I x)] + 15309 Log[1 - E^(3 I x)] + 
       109375 Log[1 - E^(-5 I x)] - 109375 Log[1 - E^(5 I x)] - 
       117649 Log[1 - E^(-7 I x)] + 117649 Log[1 - E^(7 I x)])) *)

We find four branchpoints are at

\[Beta]7 = 2 \[Pi]/{1, 3, 5, 7}
% // N

(* Out[766]= {2 \[Pi], (2 \[Pi])/3, (2 \[Pi])/5, (2 \[Pi])/7} *)

(* Out[767]= {6.28319, 2.0944, 1.25664, 0.897598} *)

The last one position is less that 1 so that a jump occurs below 1.

Now that we have understood why at n=7 something new happens and what it is, we need to explore the consequences.

The numerical consequence is the most obvious, and has been pointed out by several contributers here (see calculations above).

It remains to be made plausible that the structure of the sum s7 at x=1 is so different from that below n=7.

My plan is to invert the argument asking why the result is so simple for small n, and identify the cancelling mechanism.

Third answer

It's not a "miracle" in mathematics, it's a bug in Mathematica !

Introduction

This can be shown even simpler ! The results of the second answer remain valid, but the calculation is now much shorter, and it is left to Mathematica (with the usual level of assistance from the user).

This answer is completely sufficient to prove the existence of the bug by calculating the exact value of the sum for n=7. As the most probable cause of the bug we identify the simplification procedure.

A short list of the correct results is added for reference.

Calculations

The trick is again to introduce a parameter x writing Sin[x k] instead of Sin[k], and letting x to 1 in the end.

In fact, Mathematica immediately finds the sum with a parameter x:

s7a = Sum[(Sin[k x]/k)^7, {k, 1, \[Infinity]}]

(* Out[40]= 1/128 I (35 PolyLog[7, E^(-I x)] - 35 PolyLog[7, E^(I x)] - 
   21 PolyLog[7, E^(-3 I x)] + 21 PolyLog[7, E^(3 I x)] + 
   7 PolyLog[7, E^(-5 I x)] - 7 PolyLog[7, E^(5 I x)] - 
   PolyLog[7, E^(-7 I x)] + PolyLog[7, E^(7 I x)]) *)

In order to simplify this we use Series[] as before, giving

Series[s7, {x, 0, 10}] // Normal

(* Out[46]= (5887 \[Pi] x^6)/23040 - x^7/2 + (
 7 \[Pi] x^6 Floor[(\[Pi] - Arg[-((-1 + E^(I x))/x)] - Arg[x])/(
   2 \[Pi])])/9216 - 
...
  + (
 117649 \[Pi] x^6 Floor[(\[Pi] - 
    Arg[(E^(-7 I x) (-1 + E^(7 I x)))/x] - Arg[x])/(2 \[Pi])])/46080 *)

here we have dropped some terms from display.

Now we can simplify the expression finally with the result

Simplify[Series[s7, {x, 0, 10}] // Normal, x > 0]

(* Out[47]= (5887 \[Pi] x^6)/23040 - x^7/2 *)

which for x->1 gives us the correct result instead of the the buggy (former"miraculous") s[7].

This result supports strongly my suspicion of the OP that the error at n=7 is due to a switch in the strategy of simplification of the correct intermediate result with the PolyLogs.

It should now be easy for the Wolfram crew to detect the weak point and change it, perhaps even in the way presented here.

Summary of correct results

The results are given here in the symmetric form

$$i_s(n)=\int_{-\infty}^\infty \left(\frac{\sin (k)}{k}\right)^n \,dx$$

$$s_s(n)=\sum_{k=-\infty}^\infty \left(\frac{\sin (k)}{k}\right)^n$$

We have $i_s(n)=s_s(n)$, i.e. sum and integral are equal. They are related to the quantities used here before by $s_s(n) = 2 s(n) -1, i_s(n) = 2 i(n)$

The first 10 values are

$$\left\{\pi ,\pi ,\frac{3 \pi }{4},\frac{2 \pi }{3},\frac{115 \pi }{192},\frac{11 \pi }{20},\frac{5887 \pi }{11520},\frac{151 \pi }{315},\frac{259723 \pi }{573440},\frac{15619 \pi }{36288}\right\}$$

This series is contained in the OEIS (https://oeis.org/) as:

A002297 Numerator of (2/Pi)*Integral_{0..inf} (sin x / x)^n dx.
A002298 Denominator of (2/Pi)*Integral_{0..inf} (sin x / x)^n dx.

Here further values can be found.

Second answer

Introduction and summary

Although my first answer is just one hour ago, I need to revise myself, and come back to what my intuition was from the beginning: it is definitely a bug !

In order to show this let us just calculate the sums with a transparent method devised/used by robjohn in a slightly different context here 3

To begin with, we repeat the definitions

sd[n_] := 1 + Sum[(Sin[k]/k)^n, {k, 1, \[Infinity]}]

We use sd[n] to designate the mathematical definition while s[n] stands for the result of the Mathematica calculation.

The related integral is

i[n_] := 1/2 + Integrate[(Sin[k]/k)^n, {k, 0, \[Infinity]}]

We call sc[n] the sums calculated in the following. And we shall show that

  1. sc[n] = i[n] for n = 1..10 (and probably for any n>10)

  2. s[k] != sc[k] for n = 7 (and probably for any n>7)

To make the bug manifest the case n=7 is sufficient.

Also this disenchants the presumed "miracle": it is not a miracle but a bug in Mathematica. Most probably (and naturally) also the statement in Wolfram Mathworld (see first answer) has the same root.

The calculations

We shall consider the more general sums

s1[n_, x_] := Sum[(Sin[k x]/k)^n, {k, 1, \[Infinity]}]

They are related to the sums defined in the OP by

s[n] == 1 + s1[n, x] /. x -> 1 (* not a MMA formula *)

For n=1 there should be no doubt that

s1[1, x]

(* Out[220]= 1/2 I (Log[1 - E^(I x)] - Log[E^(-I x) (-1 + E^(I x))]) *)

I use Series[] to simplify this expression (maybe some else can find the right simplification command)

Series[%, {x, 0, 5}] // Normal

(* Out[221]= \[Pi]/2 - x/2 *)

and

sc[1] = 1 + % /. x -> 1

(* Out[222]= 1/2 + \[Pi]/2 *)

This is in agreement with s1:

s[1] // Expand

(* Out[223]= 1/2 + \[Pi]/2 *)

Generalization of the procedure to higher powers is obvious.

As an example let us take n=3

First we need the 2nd derivative of a summand

D[((Sin[k x]/k))^3, {x, 2}] // Simplify

(* Out[190]= -((3 (Sin[k x] - 3 Sin[3 k x]))/(4 k)) *)

and now execute the sum

Sum[%, {k, 1, \[Infinity]}]

(* Out[191]= 3/8 I (-Log[1 - E^(I x)] + Log[E^(-I x) (-1 + E^(I x))] + 
   3 Log[1 - E^(3 I x)] - 3 Log[E^(-3 I x) (-1 + E^(3 I x))] *)

which is in fact

h = Series[%, {x, 0, 5}] // Normal

(* Out[192]= (3 \[Pi])/4 - 3 x *)

Now we integrate twice with respect to x

Integrate[h, x]

(* Out[182]= 3/4 (\[Pi] x - 2 x^2) *)

Integrate[%, x]

(* Out[183]= 3/4 ((\[Pi] x^2)/2 - (2 x^3)/3) *)

We can simplify this double integration as a simple Mathematica construct

h = (3 \[Pi])/4 - 3 x;
Do[h = Integrate[h, x], {2}]; h

(* Out[227]= 3/4 ((\[Pi] x^2)/2 - (2 x^3)/3) *)

Finally we let x->1 to obtain

sc[3] = 1 + % /. x -> 1 // Expand

(* Out[228]= 1/2 + (3 \[Pi])/8 *)

This is identical to s3

s[3] // Expand

(* Out[229]= 1/2 + (3 \[Pi])/8 *)

i[3]

(* Out[230]= 1/2 + (3 \[Pi])/8 *)

In the same manner we retrieve the agreed results for n=4 to 6.

Now we turn to the "magic" n = 7

The 6th derivative

D[(Sin[k x]/k)^7, {x, 6}] // Simplify

(* Out[231]= (7 (-5 Sin[k x] + 2187 Sin[3 k x] - 15625 Sin[5 k x] + 
   16807 Sin[7 k x]))/(64 k) *)

The sum

Sum[%, {k, 1, \[Infinity]}]

(* Out[232]= 7/128 I (-5 Log[1 - E^(I x)] + 5 Log[E^(-I x) (-1 + E^(I x))] + 
   2187 Log[1 - E^(3 I x)] - 2187 Log[E^(-3 I x) (-1 + E^(3 I x))] - 
   15625 Log[1 - E^(5 I x)] + 15625 Log[E^(-5 I x) (-1 + E^(5 I x))] + 
   16807 Log[1 - E^(7 I x)] - 16807 Log[E^(-7 I x) (-1 + E^(7 I x))])*) 

The simplification

h = Series[%, {x, 0, 5}] // Normal

(* Out[233]= (5887 \[Pi])/32 - 2520 x *)

Now the 6-fold integration

Do[h = Integrate[h, {x, 0, y}] /. y -> x, {6}]; h

(* Out[234]= (5887 \[Pi] x^6)/23040 - x^7/2 *)

The result

sc[7] = 1 + % /. x -> 1

(* Out[235]= 1/2 + (5887 \[Pi])/23040 *)

is identical to

i[7]

(* Out[236]= 1/2 + (5887 \[Pi])/23040 *)

but clearly different from the direct Mathematica result in question

s7 = s[7]

(* Out[239]= 1 + (-23040 + 129423 \[Pi] - 201684 \[Pi]^2 + 144060 \[Pi]^3 - 
  54880 \[Pi]^4 + 11760 \[Pi]^5 - 1344 \[Pi]^6 + 64 \[Pi]^7)/46080 

as can be seen easily from the structure but finally from the numerical values

N[1/2 + (5887 \[Pi])/23040, 10] - N[s7, 10]

(* Out[241]= -9.249*10^-6 *)  

End of calculations

Contest on the magic 7

Unexpected a tacit contest was installed by two members of the community involved in this discussion, and I would like to present the first results

"Magic number 7" is a PrimeQ[7] :) – Mariusz Iwaniuk

I believe that 7 is magic, because it is the first integer greater than 2 Pi. - bbgodfrey

I would appreciate if others would join in.

First answer, 10.10.17 15:00

There has been a vivid discussion about this post in which I have tentatively asserted to have found a bug in an infinite sum. I'd like to summarize the status.

First of all: I agree now that it is not a bug, but that indeed there's a "miracle" happening at n=7.

My original conviction was based on the fact that the first 6 values of the sum had the simple structure $a + \pi b$ with rational $a$ and $b$, so why should n=7 have suddenly a much more complicated structure? and, furthermore, that these results are obtained by naively applying the Euler-Maclaurn formula, and - most of all - that I tend not to believe in miracles.

Three arguments convinced me of my error

  1. The numerical calculation of John Doty who supplemented his solution in a comment with the information that the calculation was done with a high degreee of accuracy

  2. The related discussion in MSE 1 in which the hint was given that indeed a "miracle" happens at n=7.

  3. "Proof of authority" 2: Eric Weisstein writes "Amazingly, the pattern of these sums being equal to -1/2 plus a rational multiple of pi breaks down for the power 7." and he gives the result for s[7] in my OP which Mathematica provides.

Still: Some explanation of the "Miracle" and the related "Magic number 7", some reason which should be easy to grasp would be nice to have.

References

1 https://math.stackexchange.com/questions/331404/how-to-prove-this-identity-pi-sum-limits-k-infty-infty-left-frac-sin

2 http://mathworld.wolfram.com/SincFunction.html

3 https://math.stackexchange.com/questions/453198/show-that-int-0-infty-frac-sin3xx3dx-frac3-pi8/453260#453260

$\endgroup$
  • $\begingroup$ "Magic number 7" is a PrimeQ[7] :) $\endgroup$ – Mariusz Iwaniuk Oct 10 '17 at 12:49
  • $\begingroup$ I believe that 7 is magic, because it is the first integer greater than 2 Pi. $\endgroup$ – bbgodfrey Oct 10 '17 at 13:02
  • $\begingroup$ See my answer. I think Mathematica is right here, and your Series approximation around $x=0$ and evaluation at $x=1$ is not working. $\endgroup$ – Carl Woll Oct 10 '17 at 23:49
  • $\begingroup$ @ bbgodfrey Is there someone around who can provide the results of Maple? (I don't know why I ask this only now?) $\endgroup$ – Dr. Wolfgang Hintze Oct 11 '17 at 2:55
  • $\begingroup$ @Mariusz Iwaniuk Many thanks. You have found the same "magic result" in Maple, so I surrender completely with the presumed bug, and turn to the attempt of explaining the magic7. Thinking of my "naive" usage of the Euler-Maclaurin formula which controls the difference between sum and integral I propose here another tentative "explanation": 7 is "magic" because $B_{2*7}$ is the first Bernoulli number whose absolute value is greater than 1. $\endgroup$ – Dr. Wolfgang Hintze Oct 11 '17 at 9:29
5
$\begingroup$

Using 11.1.1.

Compare with an explicitly numerical method, no symbolic analysis. To use such a method, we must define f[0]:

f[0] = 1

The sum converges pretty fast, so 100 terms is plenty:

s100[n_] := Total[Table[f[k]^n, {k, 0, 100}]]

Now try it:

Table[{n, s[n] // N, s100[n] // N}, {n, 1, 10}] // Chop

(* {{1, 2.0708, 2.06043}, {2, 2.0708, 2.06582}, {3, 1.6781, 1.6781}, {4, 1.5472, 1.5472}, {5, 1.44084, 1.44084}, {6, 1.36394, 1.36394}, {7, 1.30272, 1.30272}, {8, 1.2532, 1.2532}, {9, 1.21235, 1.21235}, {10, 1.17837, 1.17837}} *)

So, the numerical result is insensitive to the transition at n=7. I conclude that this is not a bug.

$\endgroup$
  • 1
    $\begingroup$ The results in the answer immediately above also hold true for V11.2. I would agree that there is no bug, even though the symbolic solutions become more complicated, beginning with `n == 7``. $\endgroup$ – bbgodfrey Oct 10 '17 at 0:39
  • $\begingroup$ @ John Doty I cannot follow your conclusion. You have just shown that Mathematica calculates consistently wrong values from n=7 upwards. BTW you need to take more digite to see the difference as I have shown in my OP. $\endgroup$ – Dr. Wolfgang Hintze Oct 10 '17 at 7:30
  • $\begingroup$ @bbgodfrey Did you take into account the structural Change in the symbolic expression? You can't reduce the polynomial in Pi to the structure a + Pi b with a dn b rational. I have given the proof of this staement numerically in my OP. See my comment to Akku14. $\endgroup$ – Dr. Wolfgang Hintze Oct 10 '17 at 7:38
  • $\begingroup$ @Wolfgang The errors appear to be due to routine things. s100 represents a truncated sum. As n increases, the sum converges faster, so s100 is a better approximation. At the same time, the more complicated the expression, the worse N behaves, especially in machine precision. Try Table[{n, N[s[n] - s100[n], 50]}, {n, 1, 10}]. The differences decline monotonically, as you'd expect from the faster convergence. $\endgroup$ – John Doty Oct 10 '17 at 7:44
5
$\begingroup$

This is not a bug. Mathematica calculates the sum very consequently.

The fact, that the structure of the sum simplifies for j < 7 is pure chance.

Divide the structure into parts.

(su1 = Sum[#, {k, 1, \[Infinity]}] & /@ 
         Table[TrigReduce[(Sin[k]/k)^j], {j, 1, 7, 2}] // 
    Expand) // TableForm

enter image description here

FullSimplify gives simple forms for j < 7

su1 // FullSimplify // TableForm

enter image description here

The same for even j

(su2 = Sum[#, {k, 1, \[Infinity]}] & /@ 
        Table[TrigReduce[(Sin[k]/k)^j], {j, 2, 8, 2}] // 
     Expand) // TableForm

enter image description here

su2 // FullSimplify // TableForm

enter image description here

$\endgroup$
  • $\begingroup$ @ Akku14 Thank you. This is a very useful contribution. It shows that Mathematica calculates symbolically consistent across n = 8 in terms of PolyLog functions. Hence you have narrowed down the spot of the bug: it is in the implementation of PolyLog. For n=8 the symbolic simplification of the expression gives a wrong result: a polynomial in Pi which does not reduce to the correct form a + Pi b with a and b rational: 1/2 + (151 [Pi])/630. The numerical difference appears already in the thirs decimal place (correct: 1.252984906 wrong:1.2531989845). $\endgroup$ – Dr. Wolfgang Hintze Oct 10 '17 at 7:56
  • $\begingroup$ . @Dr. Wolfgang Hintze . I don't find the symbolic simplification gives wrong results: for n=8 calculate the numerical value of the PolyLog expression before and after FullSimplify and in both cases you get 1.253198984, which you also get with NSum. I think like march mentioned, that in this case the simplified form of the Euler-Maclaurin-formula doesn't hold and the remainder R does not vanish. $\endgroup$ – Akku14 Oct 10 '17 at 8:50
  • $\begingroup$ Yes, this shows that Mathematica gives wrong values consistently, and this is due to the bug in PolyLog (symbolically as well as numerically). $\endgroup$ – Dr. Wolfgang Hintze Oct 10 '17 at 9:45
  • $\begingroup$ @Wolfgang My answer above does not use PolyLog for the truncated sum. The truncated sum is just the Total of a list of sinc functions. It's about as independent a method as you can get, but agrees to a few parts in 10^29 to the infinite Sum with n=10, truncating after 1000 terms. So this is not a bug in PolyLog. It's not a bug at all. $\endgroup$ – John Doty Oct 10 '17 at 11:06
  • $\begingroup$ @ John Ok, now that you have increased the number of valid digits appreciably and come to your previous conclusion this is convincing (assuming that Mathematica calculates correctly the powers of sinc and the additions of 1000 terms). In the meantime I have found a reference which points at an error on my side. I shall give a brief account in my answer. $\endgroup$ – Dr. Wolfgang Hintze Oct 10 '17 at 11:28
3
$\begingroup$

Here is Mathematica's answer:

m = 1 + Sum[(Sin[k]/k)^7, {k, 1, Infinity}]

1 + 1/128 I (35 PolyLog[7, E^-I] - 35 PolyLog[7, E^I] - 21 PolyLog[7, E^(-3 I)] + 21 PolyLog[7, E^(3 I)] + 7 PolyLog[7, E^(-5 I)] - 7 PolyLog[7, E^(5 I)] - PolyLog[7, E^(-7 I)] + PolyLog[7, E^(7 I)])

Here is your proposed approximation to the sum:

em = 1/2 + Integrate[(Sin[k]/k)^7, {k, 0, \[Infinity]}]

1/2 + (5887 π)/23040

Here is the difference:

N[m - em, 20]

9.2486595733354227567*10^-6 + 0.*10^-26 I

Now, as in your self-answer, consider the generalized sum:

s[x_] = 1 + Sum[(Sin[k x]/k)^7, {k, Infinity}]

1 + 1/128 I (35 PolyLog[7, E^(-I x)] - 35 PolyLog[7, E^(I x)] - 21 PolyLog[7, E^(-3 I x)] + 21 PolyLog[7, E^(3 I x)] + 7 PolyLog[7, E^(-5 I x)] - 7 PolyLog[7, E^(5 I x)] - PolyLog[7, E^(-7 I x)] + PolyLog[7, E^(7 I x)])

Notice that the PolyLog expressions come in pairs, and in particular consider the pair with $e^{-7 i x}$:

pl[x_] := I/128 (PolyLog[7, Exp[7 I x]] - PolyLog[7, Exp[-7 I x]])

Take the Series of pl, assuming $x>0$:

ser = Normal @ Series[pl[x], {x, 0, 15}, Assumptions -> x>0]

-((π^6 x)/8640) + (343 π^4 x^3)/34560 - (16807 π^2 x^5)/46080 + ( 117649 π x^6)/92160 - (117649 x^7)/92160

Take the difference of these two expressions at $x=1$:

pl[1`20] - ser /. x->1`20

9.248659573335*10^-6 + 0.*10^-22 I

Notice that this is the same as the difference that you are seeing between Mathematica's answer and your proposed approximation. So, the question is whether:

  1. The Series of pl[x] is wrong
  2. Mathematica's evaluation of pl[1`20] is wrong
  3. Something else is going on

I believe the answer is 3. Basically, I think there is some non-analyticity in the PolyLog function that Series is missing. For example, the naive MacLaurin series expansion of the non-analytic $e^{-1/x}$ at $x=0$ is just 0. Yet clearly $e^{-1/x}$ is not 0 at $x=1$. Similarly, the Series expansion of the PolyLog function has some non-analyticity that causes the series approximation at x=0 to have an incorrect value when evaluated at $x=1$.

$\endgroup$
  • $\begingroup$ @ Carl Woll I fear you are right. The non analyticity can be seen in the 6th derivative. Its value at x=1 is by a factor 11 larger than I found by expanding around x=0. See the EDIT to my third answer. $\endgroup$ – Dr. Wolfgang Hintze Oct 11 '17 at 0:46
3
$\begingroup$

This is only a copy of answer Carl Woll with the difference is that it is calculated in:

Maple 2017.3.

Here is Maple's answer:

m := 1+sum(sin(k)^7/k^7, k = 1 .. infinity)

1+(1/128*I)*(polylog(7, exp(7*I))-7*polylog(7, exp(5*I))+21*polylog(7, exp(3*I))-35*polylog(7, exp(I))+35*polylog(7, exp(-I))-21*polylog(7, exp(-3*I))+7*polylog(7, exp(-5*I))-polylog(7, exp(-7*I)))

Here is your proposed approximation to the sum:

em := 1/2+int(sin(k)^7/k^7, k = 0 .. infinity)

1/2+(5887/23040)*Pi

Here is the difference:

Digits := 30; evalf(m-em)

9,24865957333542275670*10^(-6)+1,56250000000000000000*10^(-31) I

Now, as in your self-answer, consider the generalized sum:

s := proc (x) options operator, arrow; 1+Sum(sin(x*k)^7/k^7, k = 1 .. infinity) end proc;
s1 := `assuming`([value(s(x))], [x > 0]);

1+(1/128*I)*(polylog(7, exp((7*I)*x))-7*polylog(7, exp((5*I)*x))+21*polylog(7, exp((3*I)*x))-35*polylog(7, exp(I*x))+35*polylog(7, exp(-I*x))-21*polylog(7, exp(-(3*I)*x))+7*polylog(7, exp(-(5*I)*x))-polylog(7, exp(-(7*I)*x)))

Notice that the PolyLog expressions come in pairs, and in particular consider the pair with $e^{-7 i x}$:

pl := proc (x) options operator, arrow; ((1/128)*I)*(polylog(7, exp((7*I)*x))-polylog(7, exp(-(7*I)*x))) end proc;

Take the Series of pl, assuming $x>0$:

ser := convert(series(pl(x), x = 0, 15), polynom);

-(1/8640)*Pi^6*x+(343/34560)*Pi^4*x^3-(16807/46080)*Pi^2*x^5+(1/128*I)*((117649/720)ln(-Ix)-(117649/720)ln(Ix))*x^6-(117649/92160)*x^7

Take the difference of these two expressions at $x=1$:

evalf(pl(1)-(eval(ser, x = 1)))

9,24865957333542275670*10^(-6)+0,00000000000000000000*10^0 I

and the rest(text) is the same.

EDITED:

And something from me:

with the x-dependent sum

evalf(eval(s1, x = 1))

Its value at x=1 is

1.30272435072915288380945728698+1.56250000000000000000000000000*10^(-31)*I

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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