Unexpected result with sums in set delayed functions

Question

I don't know if I've just found a bug (and if so, whether there's a workaround), or if I'm just losing my mind. Without further ado, here is my code (really, a minimal working example), run on Mathematica 12.0:

In[1]:= F[j_] := (r^2 + z^2 - 1) Sum[(-1)^k Binomial[2 j + 2, 2 k + 1] z^(2 (j - k) + 1) r^(2 k), 
    {k, 0, j}]

In[2]:= G[j_] := Sum[(-1)^(k + 1) F[k], {k, 0, j}]/(
 4 (j + 1) (j + 2) (2 j + 3))

In[3]:= G[3]

Out[3]= 1/720 (-2 z (-1 + r^2 + z^2) + (-1 + r^2 + z^2) (-4 r^2 z + 
      4 z^3) - (-1 + r^2 + z^2) (6 r^4 z - 20 r^2 z^3 + 6 z^5) + (-1 +
       r^2 + z^2) (-8 r^6 z + 56 r^4 z^3 - 56 r^2 z^5 + 8 z^7))

In[4]:= G[j]

Out[4]= -(((1 + (-1)^j + 2 r^2 - 2 (-1)^j r^2 + r^4 + (-1)^j r^4 - 
    6 r^(4 + 2 j) - 2 j r^(4 + 2 j) + 2 r^(8 + 2 j) + 
    2 j r^(8 + 2 j)) z (-1 + r^2 + z^2))/(
 4 (1 + j) (2 + j) (3 + 2 j) (-1 + r^2)^2 (1 + r^2)^2))

In[5]:= % /. j -> 3

Out[5]= -(((4 r^2 - 12 r^10 + 8 r^14) z (-1 + r^2 + z^2))/(
 720 (-1 + r^2)^2 (1 + r^2)^2))

The point of Out[5] is just to drive home the already-obvious point that Out[4] is not going to be correct. Indeed, it's not just 3: if I plug any (specific) value of j in, I get a correct answer that does not agree with G[j]. And why does G[j] have r's in the denominator anyway?!

I keep wondering if I'm just missing something obvious, which sometimes happen when I stare at code for a very long time.

Any help will be greatly appreciated.

Answer 1

To avoid unintended variable capture, one can use

F[j_]:=Module[{k},
         (r^2+z^2-1) Sum[(-1)^k Binomial[2 j+2,2 k+1] z^(2 (j-k)+1) r^(2 k),{k,0,j}]];
G[j_]:=Module[{k},
          Sum[(-1)^(k+1) F[k],{k,0,j}]/(4 (j+1) (j+2) (2 j+3))];

for example. Now

G[3]-(G[j]/.{j->3})//Simplify
(* 0 *)

---

Concerning r in the denominator, consider the much simpler case of a geometric sum:

Sum[r^k,{k,0,j}]
(* (-1+r^(1+j))/(-1+r) *)

Of course it is a polynomial for every fixed positive integer j.

Answer 2

I decided to move my comment from @user293787 to an answer as it made the comments a bit cluttered and made it unlikely that someone else encountering this issue would read it. I emphasize however that the answer that follows relies on user293787's explanation that there was conflict between the two k variables.

A simpler case where that happens :

 Clear[F,G];

 F[j_]:=Sum[j,{k,0,j}];

 G[j_]:=Sum[F[k],{k,0,j}];

 G[3] == G[j]/.j->3

output: (* False *)

TL;DR quick solution : You can use F[j_]= instead of F[j_]:= to avoid that problem.

[Edit: added the solution below that localizes the iterator in a convenient way.

More general solution : we can automatically localize all the iterators by defining a new sum and restricting the name choice for the iterators to belong to a set of one or more names.

Defining a sum that localizes the iterator k :

unique = Inactive@*Unique;
local = k -> unique[k];
sum = ReplaceAll[#, Activate@local] &@*Inactive[Sum]

One could probably include more iterators in "local" using local={k1->unique[k1],k2-> unique[k2],...}. Then instead of using Sum the user may define F and G using the new sum function :

    Clear[F,G]
F[j_]=sum[j,{k,0,j}];
G[j_]= sum[F[k],{k,0,j}];

We can check that the iterators are correct by evaluating G on a symbolic argument:

Finally we may check that the problem has been solved using sum:

G[t] == G[3] // Activate // ReplaceAll[t -> 3]

 (*True*)

End of Edit]

TL;DR Origin of the problem: symbolic handling of F[k] in G[j] where k appears in the argument of F and also in the iterator of the sum in the definition of F.

The solution above that consists of replacing := with = evades the double presence of k by simply evaluating the sum so that k no longer appears (see discussion below if difference between := and = is not clear (nothing fancy)).

The symbolic handling occurs when the upper bound for the iterator in G (the last element of the second argument in Sum) is symbolic but not when it is an integer. This leads to a difference between G[3] and G[j]/.j->3. The solution

The problem is tricky and interesting and it could be quite difficult to debug in a long code.

Explanation

Pictorial summary and reference for the discussion below using traceView2 from this answer (the images and discussion below will focus on the definitions of F and G of the simpler example above and should not be confused with the original F and G from the original question. Also, the discussion below does not require the reader to understand the images below):

Consider first G[3].

From traceview2[G[3]] (output shown in first picture), one sees that the evaluation process is as expected:

G[3]=F[0]+F[1]+F[2]+F[3]

Consider now G[j] with a symbolic j. Mathematica does not know the upper bound in this case so it can not do:

F[0]+F[1]+...

Instead, as can be checked using traceView2[G[j]] (second image), Mathematica attempts to do pure symbolic manipulations replacing F[k] by the symbolic expression it gets from the definition of F, that is, the output of F[k], but...

F[k]

evaluates to

Sum[k,{k,0,k}]

which according to Mathematica is k(1+k)/2 but according to Mathematics, it is an invalid request.

The issue there is that maybe Mathematica applies the formula Sum[k,{k,0,n}] then takes n=k at the end where k is a global variable.

So,

Computing F[j], we obtain the expected result j(j+1) and then replacing j by k with /.j->k, we find a result that is different than the result obtained from F[k].

---

Ok, but why did that not happen with G[3] ?

Attributes[Sum]

output: {HoldAll, Protected, ReadProtected}

For the following mini-discussion that answers the above question,the important component in the above output is HoldAll.

Sum holds F[k] meaning F[k] is not replaced by its definition right away. Instead at the beginning, basically, Mathematica sees the string of characters "F" "[" "k" "]" without checking what is under the hood (at least conceptually not sure what Mathematica is actually doing). Then using the iterator k in the second argument of Sum, Sum successively makes replacements for k, k->0; F[0], k->1; F[1], k->2; F[2], k->3; F[3] and each F[0],F[1],F[2],F[3] will be replaced by its definition (I do not know if replacement rules are actually being used).

When the upper bound is symbolic that sequence can not happen and so F[k] is evaluated with a symbolic k

---

Solution:

If instead F[j_]=Sum[k,{k,0,j}] (no :=) then there are no surprises because the sum is evaluated before assigning the expression to the function F. Hence, Mathematica sees F[j_]=j(j+1).

So a possible solution for the initial question of OP would be (replaced F[j_]:= with F[j_]=) :

Clear[F,G];
F[j_]=(r^2+z^2-1) Sum[(-1)^k Binomial[2 j+2,2 k+1] z^(2 (j-k)+1) r^(2 k),{k,0,j}];
G[j_]:=Sum[(-1)^(k+1) F[k],{k,0,j}]/(4 (j+1) (j+2) (2 j+3));

Appendix: Difference between F[j_]= and F[j_]:= and when to use one or the other (nothing fancy in the discussion below).

If I remember well, for a somewhat long time working with Mathematica I thought I always had to define a function with the construct F[j_]:=expression.

However, often it is better to use F[j_]= either to avoid surprises like in the above example or to increase the speed of execution as F[j_]:= means nothing is being evaluated until it is called with an argument. For example, compare (* The second one takes 2 seconds before it outputs something*):

Clear[f];AbsoluteTiming[f[j_]:=(Pause[2];Print["that took a long time"])]

Clear[f];AbsoluteTiming[f[j_]=(Pause[2];Print["that took a long time"])]

I use f[j_]:= when a symbolic expression can not be obtained or I need a black box numerical function for certain numerical computations in which case I would use f[j_?NumericQ]:=.