A case when the sum rule in integration does not work in Mathematica

Question

Bug introduced in 8.0 or earlier and persisting through 11.2.0 or later. Fixed in 13.2.0 or earlier.

---

I am using Mathematica 8.0.1 and discovered that this integral:

a = Sum[Integrate[1/n^s*2^(s - 1), s], {n, 1, 6}]

$$\frac{s}{2}-\frac{2^{-s-1}}{\log (2)}+\frac{2^{s-1}}{\log (2)}+\frac{2^{s-1} 3^{-s}}{\log (2)-\log (3)}-\frac{3^{-s}}{2 \log (3)}+\frac{2^{s-1} 5^{-s}}{\log (2)-\log (5)}$$

and this integral:

b = Integrate[Sum[1/n^s*2^(s - 1), {n, 1, 6}], s]

$$\frac{1}{2} \left(\frac{\left(\frac{2}{3}\right)^s}{\log (48)}-\frac{2^{-s}}{\log (2)}+\frac{2^s}{\log (144)}+\frac{3^{-s}}{\log \left(\frac{10}{3}\right)}+\frac{\left(\frac{2}{5}\right)^s}{\log \left(\frac{288}{5}\right)}\right)$$

give symbolically different results although they according to the sum rule in integration should be the same.

The sum rule has been discussed earlier on this forum but I could not find a question that would answer this particular example.

The same difference is seen in this numeric test:

(*Start*)
 Clear[a, b, s, n];
 a = Sum[Integrate[1/n^s*2^(s - 1), s], {n, 1, 6}]
 b = Integrate[Sum[1/n^s*2^(s - 1), {n, 1, 6}], s]
 s = 2;
 Print["The following difference should be zero:"]
 N[a - b]
 (*End*)

where the numeric value 2.67373 from the output, should be 0.

Is this a bug / unexpected behaviour in other versions of Mathematica too?

I have been in contact with the Technical Support at Wolfram Research Inc., and they agree that this difference does seem incorrect.

Edit this question if you like.

---

Edit 11.2.2018:

Compare this integral:

a = Integrate[1/1^s*2^(s - 1), s] + Integrate[1/2^s*2^(s - 1), s] + Integrate[1/3^s*2^(s - 1), s]

$$\frac{s}{2}+\frac{2^{s-1}}{\log (2)}+\frac{2^{s-1} 3^{-s}}{\log (2)-\log (3)}$$

to this integral:

b = Integrate[1/1^s*2^(s - 1) + 1/2^s*2^(s - 1) + 1/3^s*2^(s - 1), s]

$$\frac{1}{2} \left(s-\frac{\left(\frac{2}{3}\right)^s}{\log \left(\frac{3}{2}\right)}+\frac{2^s}{\log (4)}\right)$$

The latter is wrong.

Answer 1

This is a report that the are no issues in the latest -13.2.0- version.

With a clen kernel

a = Sum[Integrate[1/n^s*2^(s - 1), s], {n, 1, 6}];
b = Integrate[Sum[1/n^s*2^(s - 1), {n, 1, 6}], s];
a - b // FullSimplify

Also,

Clear[a, b, s, n];
a = Sum[Integrate[1/n^s*2^(s - 1), s], {n, 1, 6}];
b = Integrate[Sum[1/n^s*2^(s - 1), {n, 1, 6}], s];
s = 2;
Print["The following difference should be zero:"];
N[a - b]

returns

And finally,

a = Integrate[1/1^s*2^(s - 1), s] + Integrate[1/2^s*2^(s - 1), s] + 
   Integrate[1/3^s*2^(s - 1), s];
b = Integrate[1/1^s*2^(s - 1) + 1/2^s*2^(s - 1) + 1/3^s*2^(s - 1), 
   s];
a - b // FullSimplify

For clarity and completeness, I am attaching a screenshot

Answer 2

Instead of comparing with the sum rule, it is perhaps more direct to differentiate the integral and compare. When integrating just two terms of the sum for n up to 10, there are 7 pairs which integrate incorrectly:

ClearAll[test];
mem : test[i_, j_] := mem = Module[{func, prim},
    func = 1/i^s*2^(s - 1) + 1/j^s*2^(s - 1);
    prim = Integrate[func, s];
    Simplify[func - D[prim, s]] === 0
    ];

nn = 10;
failed = Flatten@Table[! test[i, j], {i, nn}, {j, i + 1, nn}];
Pick[Flatten[Table[{i, j}, {i, nn}, {j, i + 1, nn}], 1], failed]

(*  {{1, 6}, {1, 10}, {3, 5}, {3, 7}, {5, 7}, {5, 9}, {7, 9}}  *)

Checking the first failure explicitly:

With[{i = 1, j = 6},
  func = 1/i^s*2^(s - 1) + 1/j^s*2^(s - 1)];
prim = Integrate[func, s]
func - D[prim, s] // Simplify
func - D[prim, s] /. s -> 1.
(*
  1/2 (-(3^-s/Log[3]) + 2^s/Log[12])

  (2^(-1 + s) Log[6])/Log[12]
  0.721057
*)

Answer 3

These integration results are the same in version 11.2.

It seems as though the issue may be related the indefinite integration. The expressions produced if an arbitrary definite integration is chosen instead are equivalent:

a = Sum[Integrate[1/n^s*2^(s - 1), {s, si, sf}], {n, 1, 6}];
b = Integrate[Sum[1/n^s*2^(s - 1), {n, 1, 6}], {s, si, sf}];
Simplify[a==b]

True

Of course, any sum of different constants of integration would be expected to also be constant, which is not the case here.

As an even more minimal case of issues:

a = Sum[Integrate[1/n^s*2^(s - 1), s], {n, {3, 5}}];
b = Integrate[Sum[1/n^s*2^(s - 1), {n, {3, 5}}], s];
Simplify[a == b]

Produces a variable expression as well, despite only being the sum of two terms.