6
$\begingroup$

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

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?

Issue persists 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.

$\endgroup$
  • 1
    $\begingroup$ 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. :) $\endgroup$ – Michael E2 Oct 15 '15 at 22:56
  • $\begingroup$ @MichaelE2 I tried Trace but did not learn much. Thanks for the suggestion, though. $\endgroup$ – bbgodfrey Oct 16 '15 at 1:48
  • $\begingroup$ If I followed the computation correctly, Integrate mistakenly concludes that all the roots of the denominator lie on the real line: i.stack.imgur.com/peMwW.png -- no time to investigate further right now $\endgroup$ – Michael E2 Oct 16 '15 at 12:56
  • $\begingroup$ @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. $\endgroup$ – bbgodfrey Oct 16 '15 at 15:20
  • 2
    $\begingroup$ Internal`Integrate`debugSwitch. Almost 10MB of output, so beware. (My institution is slow about making available new versions, so no 10.3 for me yet.) $\endgroup$ – Michael E2 Oct 16 '15 at 17:47
4
$\begingroup$

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 - ω));

Compact answer

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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