Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
@belisarius I updated it. –  popffabrik Feb 17 at 22:18
1  
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 ... –  Szabolcs Feb 17 at 22:19
    
@Szabolcs Thanks, noted. –  popffabrik Feb 17 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! –  Szabolcs Feb 17 at 22:29
1  
Re-added the bugs tag in view of my answer. Please let me know if someone thinks otherwise –  belisarius Feb 18 at 2:32

2 Answers 2

up vote 2 down vote accepted

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}];
(* instead of 
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
*)
share|improve this answer
    
As a side note, each of the Nx gave a Complex result when evaluated with the original limits sequence –  belisarius Feb 18 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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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