# Inconsistent results for equivalent converging symbolic integrals

Bug introduced in 9.0 or earlier and persisting through 11.0.1 or later

I have looked at previous questions and I'm aware that this seems to be a known bug: Mathematica giving inconsistent results for symbolic integrals done in different ways. The origins for the bugs seem to be varied.

I'm having a hard time tracking down where Mathematica is going wrong with this particular one, especially since I am explicitly assuming that all my variables ($L$ and $\rho$) are positive. I've played around with the assumptions without any luck. Here's my code:

$Assumptions = {ρ > 0, L > 0}; limits = Sequence[{ux, 0, L}, {vx, 0, L}, {uy, 0, ux}, {vy, 0, vx}] N1 = Integrate[a1 = Exp[-ρ (ux - uy) (vx - vy)], limits]; N2 = Integrate[a2 = Exp[-ρ (ux - uy) (vx - vy)] (ρ (ux - uy) (vx - vy)), limits]; N3 = Integrate[a3 = Exp[-ρ (ux - uy) (vx - vy)] (ρ (ux - uy) (vx - vy))^2, limits]; N4 = Integrate[a4 = Exp[-ρ (ux - uy) (vx - vy)] (ρ (ux - uy) (vx - vy))^3, limits]; int1 = Integrate[4 ρ (a1 - 7/2 a2 + 2 a3 - 1/4 a4), limits] int2 = Simplify[4 ρ (N1 - 7/2 N2 + 2 N3 - 1/4 N4)] // Expand  Apologies for the large expressions but the particular simplifcations only happen when all four integrands are included with the correct factors. By comparing integrands it's easy to see that int1 and int2 are the same integrals. But for int1 I get (1 - E^(-L^2 ρ)) L^2  whereas for int2 I get L^2 - E^(-L^2 ρ) L^2 - (7 L^4 ρ)/4 + 7/2 L^4 ρ Log[L]  I know that the outcome for int1 is correct whereas that for int2 is incorrect. Any thoughts would be very welcome. • @belisarius I updated it. Feb 17 '14 at 22:18 • Please don't use the bugs tag when you post a new question. This tag has a special use and it's only added after the community (or WRI support) has agreed that it's really a bug ... Feb 17 '14 at 22:19 • @Szabolcs Thanks, noted. Feb 17 '14 at 22:20 • Actually it wasn't so easy for me to see that the two integrals are really the same (and also check that you didn't make a typo) because the code is rather large. So I edited your post a bit to make this very clear. Please check the edit! Feb 17 '14 at 22:29 • Re-added the bugs tag in view of my answer. Please let me know if someone thinks otherwise Feb 18 '14 at 2:32 ## 2 Answers It's a bug, of course. Mathematica gets dizzy by the "hanging" integration order (I believe). As simple as it is, it was tough to find out but it gives the right result if you just change the order of limits: $Assumptions = {ρ > 0, L > 0};

limits = Sequence[{ux, 0, L}, {uy, 0, ux}, {vx, 0, L}, {vy, 0, vx}];
limits = Sequence[{ux, 0, L}, {vx, 0, L}, {uy, 0, ux}, {vy, 0, vx}]
*)

N1 = Integrate[a1 = Exp[-ρ (ux - uy) (vx - vy)],                            limits];
N2 = Integrate[a2 = Exp[-ρ (ux - uy) (vx - vy)] (ρ (ux - uy) (vx - vy)),    limits];
N3 = Integrate[a3 = Exp[-ρ (ux - uy) (vx - vy)] (ρ (ux - uy) (vx - vy))^2,  limits];
N4 = Integrate[a4 = Exp[-ρ (ux - uy) (vx - vy)] (ρ (ux - uy) (vx - vy))^3,  limits];

int1 = Integrate[4 ρ (a1 - 7/2 a2 + 2 a3 - 1/4 a4), limits]

int2 = FullSimplify[4 ρ (N1 - 7/2 N2 + 2 N3 - 1/4 N4)]

(*
(1 - E^(-L^2 ρ)) L^2
(1 - E^(-L^2 ρ)) L^2
*)

• As a side note, each of the Nx gave a Complex result when evaluated with the original limits sequence Feb 18 '14 at 2:30

I have isolated a simplified instance of the bug:

a = Exp[-(ux - uy) (vx - vy)] ((ux - uy) (vx - vy))^2;
i1 = Integrate[a, {vy, 0, vx}];
Assuming[ux > 0,
i2a = Integrate[i1, {vx, 0, 1}, {uy, 0, ux}];
i2b = Integrate[Integrate[i1, {uy, 0, ux}], {vx, 0, 1}];
i2b - i2a // FullSimplify
]


which gives 4 I \[Pi] (incorrect) rather than 0 (correct).

We can trace the likely source of this bug by inspecting the first integral i1 which is (2 - Gamma[3, (ux - uy) vx])/(ux - uy) - i.e. a pair of terms each of which has a pole at uy = ux, but which have oppositely-signed contributions so they cancel out.

However, this cancellation does not occur when the subsequent integrals are evaluated as a single Integrate - i.e. i2a - where it seems that the oppositely-signed contributions become same-signed contributions.

So the bug seems to be inconsistent handling of the 2 contributions to the (illusory) pole at uy = ux when the integral is evaluated as a single Integrate.