Double integral yields zero unexpectedly

Question

In developing a large set of functions and definitions, I sometimes got zero unexpectedly. In trying to identify the problem via Trace, the fact that Integrate doesn't give much insight in its inner workings turned out unhelpful.

Eventually I realized that the issue can be formulated as a minimal working example as follows (neither x nor f[x] are defined in the below; so this is purely symbolic):

Integrate[f[x] x^2, x] - Integrate[f[x], x] x^2
Integrate[Integrate[f[x] x^2, x] - Integrate[f[x], x] x^2, x]
Integrate[Integrate[f[x] x^2, x], x] - Integrate[Integrate[f[x], x] x^2, x]

(*
Out[1]= Integrate[f[x] x^2, x] - Integrate[f[x], x] x^2
Out[2]= 0
Out[3]= Integrate[Integrate[f[x] x^2, x], x] - Integrate[Integrate[f[x], x] x^2, x]
*)

So, the first line is not zero of itself. Using it as integrand in another integration over x, however, does yield zero. In fact, it does so when replacing x^2 with x^n for all n>1.

From Integration by Parts by hand, however, I reckon that the correct answer is

-n Integrate[Integrate[f[x],x] x^(n-1), x]

integrated once more over x, which is non-zero in general. Moreover, given the linearity of Integration, I was surprised that the second and third input line do not give rise to the same result, according to Mathematica.

I suspect that source of confusion may be the multiple use of x. Indeed, for specified functions, or when the Integrals are given appropriate limits, the issue disappears. As I am interested in the Indefinite integrals, this does not solve my problem. Am I missing something, mathematically? Or is this explainable behaviour of Mathematica? And, is there a workaround, that would always give me the result Out[3]?

Answer 1

From more of a programming point of view, there is an ambiguity in

Integrate[Integrate[f[x] x^2, x] - Integrate[f[x], x] x^2, x]

The problem is of scoping. There are in essence, three variables x, each belonging to one integral. In the simpler expression

Integrate[f[x] x^2, x] - Integrate[f[x], x] x^2

it is clear that the x in each integral is distinct and bound inside the integral and that the x^2 outside the second integral is free. But wrap this in another integral that binds x and the question becomes which x's belong to the third integral.

In some sense the answer is all of them, because the results of the interior integrals are depend on x which is to be integrated by the outer integral. When properly interpreted, the interior integrals are to be evaluated first and then integrated by the outer integral. But in this case, the inner Integrate[] codes return unevaluated. It may be (and it seems likely to me) that then all the instances of x in the integrand are treated as belonging to the outer integral. It's possible that x^2 Integrate[f[x], x] and Integrate[x^2 f[x], x] are mistakenly treated as equivalent. I'm not sure how to verify this without wading through a lot of internal trace output, which might be inconclusive. In some sense, it's a bug, but it might be a difficult one to fix.

One workaround it to use different symbols for different contexts, but that has drawbacks, too:

Integrate[
 Integrate[f[x] x^2, {x, 0, y}] - Integrate[f[x], {x, 0, y}] y^2,
 {y, 0, z}]

Integrate[Integrate[f[x] x^2, {x, 0, y}], {y, 0, z}] - 
 Integrate[Integrate[f[x], {x, 0, y}] y^2, {y, 0, z}]

---

Update: It appears from an internal trace that in the integral

Integrate[Integrate[f[x] x^2, x] - Integrate[f[x], x] x^2, x]

the constant x^2 multiplying the second integral of the integrand is replaced with a placeholder via a rule x^2 -> Integrate`Y. But that also replaces the x^2 in the first integral, which indeed would result in 0.