# Unexpected Behavior from Integrate with PrincipalValue -> True

Bug introduced in 9.0 or earlier and persisting through 12.1 or later. Reported to Wolfram, Inc. as CASE:4544446.

The following code, taken from my answer to 97024,

Integrate[1/((1 - 4 ωp^2 - 4 I ωp γ) (1 - ωp^2 - 2 I ωp γ)^2 (ωp - ω)),
{ωp, -Infinity, Infinity}, Assumptions -> γ > 0 && ω > 0, PrincipalValue -> True]


runs for 41 minutes on Mathematica

$Version (* 10.2.0 for Microsoft Windows (64-bit) (July 7, 2015) *)  before returning unevaluated, in the process twice giving the two error messages, Limit::cas: Warning: contradictory assumption(s) False&&ωp>4096 encountered. >> Integrate::idiv: Integral of ... does not converge on does not converge on {-∞,∞}. >>  neither of which is valid. Omitting the term, (1 - ωp^2 - 2 I ωp γ)^2 or setting ω -> 0 produces the same outcome, only faster. In fact, the only variant I have been able to make work is to set γ to a numerical value, here 3. Integrate[1/((1 - 4 ωp^2 - 4 I ωp 3) (1 - ωp^2 - 2 I ωp 3)^2 (ωp - ω)), {ωp, -Infinity, Infinity}, Assumptions -> ω > 0, PrincipalValue -> True] (* -((I π)/((-1 + 4 ω (3 I + ω)) (-1 + ω (6 I + ω))^2)) *)  Is there some work-around for this problem, not involving evaluation of residues? Issues persist in 10.3 With $Version
(* 10.3.0 for Microsoft Windows (64-bit) (October 9, 2015) *)


a simpler version of the code above,

Integrate[1/((1 - 4 ωp^2 - 4 I ωp γ) (ωp - ω)), {ωp, -Infinity, Infinity},
Assumptions -> γ > 0 && ω > 0, PrincipalValue -> True]


produces the same spurious error messages as before, although it takes twice as long to do so (about eight minutes). Incidentally, a still faster way to reproduce this problem is

Integrate[1/(-(ωp - 1/2 (-I γ - Sqrt[1 - γ^2])) (ωp - ω)), {ωp, -∞, ∞},
Assumptions -> γ > 0 && ω ∈ Reals, PrincipalValue -> True]


which takes only seconds. Note that I also replaced ω > 0 by ω ∈ Reals to see whether it mattered. It does not.

Issues persist for 12.1

Mathematica 12.1 now can solve the first integral given under "Issues persist in 10.3" (in just over two minutes), but the answer is wrong.

(* ConditionalExpression[(Pi*(I*γ - Sqrt[1 - γ^2] + 2*ω))/(Sqrt[1 - γ^2]*
(I + 4*(γ - I*ω)*ω)), γ < 1] *)


as can be seen from summing the residues,

Residue[1/((1 - 4 ωp^2 - 4 I ωp γ) (ωp - ω)), {ωp, #}] & /@
(ωp /. Solve[((1 - 4 ωp^2 - 4 I ωp γ) (ωp - ω)) == 0, ωp])
(* {-(1/(2 Sqrt[1 - γ^2] (I γ + Sqrt[1 - γ^2] + 2 ω))),
-(1/(2 Sqrt[1 - γ^2] (-I γ + Sqrt[1 - γ^2] - 2 ω))),
I/(I + 4 γ ω - 4 I ω^2)} *)

Simplify[2 π I (Total[%] - Last[%]/2)]
(* π/(I + 4 γ ω - 4 I ω^2) *)


The other integrals given in this question either return unevaluated with error messages or run indefinitely.

• I believe the ωp>4096 comes from an internal Limit making an assumption, if it's any help Trace-ing down the problem...I'll probably delete this, as it probably won't prove to be much help. But it's all the time I got right now. :) Commented Oct 15, 2015 at 22:56
• @MichaelE2 I tried Trace but did not learn much. Thanks for the suggestion, though. Commented Oct 16, 2015 at 1:48
• If I followed the computation correctly, Integrate mistakenly concludes that all the roots of the denominator lie on the real line: i.sstatic.net/peMwW.png -- no time to investigate further right now Commented Oct 16, 2015 at 12:56
• @MichaelE2 Thanks for continuing to investigate this question. How did you obtain the PV limits? I did not obtain them with Trace. You may be correct that Integrate thinks they are on the real axis, but that would be really strange, because the expressions contain I. By the way, I installed 10.3 this morning and shall rerun this calculation with it. Quite slow, though. Commented Oct 16, 2015 at 15:20
• InternalIntegratedebugSwitch. Almost 10MB of output, so beware. (My institution is slow about making available new versions, so no 10.3 for me yet.) Commented Oct 16, 2015 at 17:47

Interesting exercise. It can be done with some preparatory steps in a few seconds of computing time.

The OP asks for the integral

q := Integrate[h, {ωp, -∞, ∞},
Assumptions -> {γ > 0, ω > 0}, PrincipalValue -> True]


where the integrand is

h = 1/((1 - 4 ωp^2 -
4 I ωp γ) (1 - ωp^2 -
2 I ωp γ)^2 (ωp - ω));


Letting

g = 1/((x - a) (x - b) (x - c) (x - d)^2 (x - f)^2);
rep = Thread[{b, c, d, d, f, f, a} -> (ωp /. Solve[0 == 1/h, ωp])];


we have

Simplify[h == (-1/4) g /. rep /. x -> ωp]

(*
Out[145]= True
*)


and the integral becomes

q = Timing[(-1/4) Simplify[
Plus @@ Limit[
Integrate[List @@ Apart[g], {x, -z, z}, PrincipalValue -> True,
Assumptions -> {z > 0, a ∈ Reals, Im[b] != 0, Im[c] != 0,
Im[d] != 0, Im[f] != 0}], z -> ∞] /. rep]]
{14.2116911, -((
2 π)/((I + 4 γ ω - 4 I ω^2) (2 γ ω -
I (-1 + ω^2))^2))}


and it is done in less that 15 seconds (Version 10.1. on a normal Windows Laptop).

Detailed derivation

We show here that the integral can in fact be calculated refraining from the use of evaluating residues. The main idea is to understand the integral as a limit z->∞ of the finite integral between -z and z. Also, an explicit partial fraction decomposition using Apart[] helps reduce the computing time.

Determining explcitly the roots of the denominator

sol = Solve[0 == 1/h, ωp]

(*
Out[70]= {{ωp -> 1/2 (-I γ - Sqrt[1 - γ^2])}, {ωp ->
1/2 (-I γ + Sqrt[1 - γ^2])}, {ωp -> -I γ -
I Sqrt[-1 + γ^2]}, {ωp -> -I γ -
I Sqrt[-1 + γ^2]}, {ωp -> -I γ +
I Sqrt[-1 + γ^2]}, {ωp -> -I γ +
I Sqrt[-1 + γ^2]}, {ωp -> ω}}
*)


Letting

p = ωp /. sol

(*
Out[88]= {1/2 (-I γ - Sqrt[1 - γ^2]),
1/2 (-I γ + Sqrt[1 - γ^2]), -I γ -
I Sqrt[-1 + γ^2], -I γ -
I Sqrt[-1 + γ^2], -I γ +
I Sqrt[-1 + γ^2], -I γ + I Sqrt[-1 + γ^2], ω}
*)


We find

Simplify[4 Times @@ (ωp - p) == -1/h]

(*
Out[89]= True
*)


Assuming general parameters defined in this replacement

rep = Thread[{b, c, d, d, f, f, a} -> p]

(*
Out[90]= {b -> 1/2 (-I γ - Sqrt[1 - γ^2]),
c -> 1/2 (-I γ + Sqrt[1 - γ^2]),
d -> -I γ - I Sqrt[-1 + γ^2],
d -> -I γ - I Sqrt[-1 + γ^2],
f -> -I γ + I Sqrt[-1 + γ^2],
f -> -I γ + I Sqrt[-1 + γ^2], a -> ω}
*)


We shall consider the integral of the more general expression for the integrand (for ease of reading we let ωp -> x)

g = 1/((x - a) (x - b) (x - c) (x - d)^2 (x - f)^2);

g /. rep /. x -> ωp // Simplify

(*
Out[97]= -(4/((ω - ωp) (-1 +
2 I γ ωp + ωp^2)^2 (-1 + 4 I γ ωp +
4 ωp^2)))
*)


The integral over g is then -4 times the integral in question, i.e. q = -g/4.

Before taking the integral It helps decomposing g into partial fractions, and splitting it into a list

ga = List @@ Apart[g]

(*
Out[107]= {1/((a - b) (a - c) (a - d)^2 (a - f)^2 (-a + x)), -(
1/((a - b) (b - c) (b - d)^2 (b - f)^2 (-b + x))), -(
1/((a - c) (-b + c) (c - d)^2 (c - f)^2 (-c + x))), -(
1/((a - d) (-b + d) (-c + d) (d - f)^2 (-d + x)^2)), (
2 a b c - 3 a b d - 3 a c d - 3 b c d + 4 a d^2 + 4 b d^2 + 4 c d^2 - 5 d^3 +
a b f + a c f + b c f - 2 a d f - 2 b d f - 2 c d f +
3 d^2 f)/((a - d)^2 (-b + d)^2 (-c + d)^2 (d - f)^3 (-d + x)), -(
1/((b - f) (-a + f) (-c + f) (-d + f)^2 (-f + x)^2)), (
2 a b c + a b d + a c d + b c d - 3 a b f - 3 a c f - 3 b c f - 2 a d f -
2 b d f - 2 c d f + 4 a f^2 + 4 b f^2 + 4 c f^2 + 3 d f^2 -
5 f^3)/((a - f)^2 (-b + f)^2 (-c + f)^2 (-d + f)^3 (-f + x))}
*)


Now the integral over the list ga is quickly calculated

gai = Timing[
Integrate[ga, {x, -z, z}, PrincipalValue -> True,
Assumptions -> {z > 0, a ∈ Reals, Im[b] != 0, Im[c] != 0,
Im[d] != 0, Im[f] != 0}]]

(*
Out[110]= {15.8029, {(-Log[-a - z] +
Log[-a + z])/((a - b) (a - c) (a - d)^2 (a - f)^2), (-Log[b - z] +
Log[b + z])/((a - b) (b - c) (b - d)^2 (b - f)^2), (-Log[c - z] +
Log[c + z])/((a - c) (-b + c) (c - d)^2 (c - f)^2), -((
2 z)/((a - d) (-b + d) (-c + d) (d - f)^2 (d^2 -
z^2))), ((d (-3 b c + 4 (b + c) d - 5 d^2) + (b c - 2 (b + c) d +
3 d^2) f +
a (2 b c - 3 b d - 3 c d + 4 d^2 + (b + c - 2 d) f)) (Log[d - z] -
Log[d + z]))/((a - d)^2 (b - d)^2 (c - d)^2 (d - f)^3), -((
2 z)/((b - f) (d - f)^2 (-a + f) (-c + f) (f^2 - z^2))), ((2 a b c +
a b d + a c d +
b c d - (3 b c + 3 a (b + c) + 2 (a + b + c) d) f + (4 (a + b + c) +
3 d) f^2 - 5 f^3) (Log[f - z] - Log[f + z]))/((a - f)^2 (b - f)^2 (c -
f)^2 (-d + f)^3)}}
*)


Now taking the limit z->∞

gail = Timing[Limit[gi[[2]], z -> ∞]]

(*
Out[111]= {0.265202, {-((I π)/((a - b) (a - c) (a - d)^2 (a - f)^2)), -((
I π)/((a - b) (b - c) (b - d)^2 (b - f)^2)), -((
I π)/((a - c) (-b + c) (c - d)^2 (c - f)^2)), 0, (
I (d (-3 b c + 4 (b + c) d - 5 d^2) + (b c - 2 (b + c) d + 3 d^2) f +
a (2 b c - 3 b d - 3 c d + 4 d^2 + (b + c - 2 d) f)) π)/((a -
d)^2 (b - d)^2 (c - d)^2 (d - f)^3), 0, (
I (2 a b c + a b d + a c d +
b c d - (3 b c + 3 a (b + c) + 2 (a + b + c) d) f + (4 (a + b + c) +
3 d) f^2 - 5 f^3) π)/((a - f)^2 (b - f)^2 (c - f)^2 (-d + f)^3)}}
*)


and summing up the terms

gs = Plus @@ gail[[2]]

(*
Out[117]= -((I π)/((a - b) (a - c) (a - d)^2 (a - f)^2)) - (
I π)/((a - b) (b - c) (b - d)^2 (b - f)^2) - (
I π)/((a - c) (-b + c) (c - d)^2 (c - f)^2) + (
I (2 a b c + a b d + a c d +
b c d - (3 b c + 3 a (b + c) + 2 (a + b + c) d) f + (4 (a + b + c) +
3 d) f^2 - 5 f^3) π)/((a - f)^2 (b - f)^2 (c - f)^2 (-d + f)^3) + (
I (d (-3 b c + 4 (b + c) d - 5 d^2) + (b c - 2 (b + c) d + 3 d^2) f +
a (2 b c - 3 b d - 3 c d + 4 d^2 + (b + c - 2 d) f)) π)/((a -
d)^2 (b - d)^2 (c - d)^2 (d - f)^3)
*)


Replacing the parameters

gsrep = gs /. rep // Simplify

(*
Out[120]= (8 π)/((I + 4 γ ω - 4 I ω^2) (2 γ ω -
I (-1 + ω^2))^2)
*)


Hence the original integral is simply

q = -1/4 gsrep
`

$$q = -\frac{2 \pi }{\left(4 \gamma \omega -4 i \omega ^2+i\right) \left(2 \gamma \omega -i \left(\omega ^2-1\right)\right)^2}$$

Done.