8
$\begingroup$

Bug introduced in 12.2 or earlier and persisting through 13.0.1.

How to calculate the following symbolic integral in MMA?

Every variable is real. Although the integrand is always positive, MMA (mine is V12.2) just gives zero as far as I've tried.
It is wrong for sure. Note also that if we set b=0 in the first place, it can give the corresponding correct finite result.

f = 1/((-(a - I b) + (x + I c)^2) (-(a + I b) + (x - I c)^2));
Integrate[f, {x, -∞, ∞}]
Integrate[f, {x, -∞, ∞}, Assumptions -> c > 0 && a > 0 && b ∈ Reals]

Edit

Taking into account several comments under this question (and anwers as well) it is plausible rewriting the integral to another form where the system returns manifestly wrong results justifying the bugsx tag. Namely the integral can be recast to the form positive constant times non-zero integral, e.g.

Integrate[ 1/(4 p^2 (q + z)^2 + (-1 - p^2 + z^2)^2), {z, -∞, ∞}, 
           Assumptions -> p > 0 && q == Sqrt[p^2 + 1]]

0

while setting for p a specific value, e.g. p = 1;

Integrate[ 1/(4 p^2 (q + z)^2 + (-1 - p^2 + z^2)^2), {z, -∞, ∞}, 
           Assumptions -> q == Sqrt[p^2 + 1]]

the system returns correctly a warning that the integral does not converge.

Improve this question
$\endgroup$
9
  • $\begingroup$ Mathematica is right in view of r = Simplify[f, Assumptions -> {a, b, c} \[Element] Reals && x > 0];s=Simplify[f,Assumptions->Assumptions->{a,b,c}\[Element]Reals&&x>0&&x<0];r-s which results in 0. $\endgroup$
    – user64494
    Commented Mar 26, 2022 at 11:28
  • $\begingroup$ @user64494 The integrand is positive.... $\endgroup$
    – xiaohuamao
    Commented Mar 26, 2022 at 11:43
  • $\begingroup$ xiaohuamao (@ does not work): The integrand is complex-valued, not positive. Master your math. $\endgroup$
    – user64494
    Commented Mar 26, 2022 at 12:10
  • $\begingroup$ @user64494 The integrand is in the form $a\times a^*$, which is positive. $\endgroup$
    – xiaohuamao
    Commented Mar 26, 2022 at 12:16
  • $\begingroup$ xiaohuamao (@ does not work): You are right, Mathematica too. $\endgroup$
    – user64494
    Commented Mar 26, 2022 at 12:32

3 Answers 3

Reset to default
7
$\begingroup$

Split the integration range from -Infinity to 0 and from 0 to Infinity. Add the two Rootsums and convert the result to radicals, then help Mathematica simplify the result. You will find

NIntegrate[1/(a^2 + b^2 + 4*b*c*x + 2*a*(c^2 - x^2) + (c^2 + x^2)^2),
{x, -100, 100}, WorkingPrecision -> 100]

equals

(2*Pi*c)/(4*c^2*(a + c^2) - b^2)

if the conditions

c > Sqrt[Sqrt[a^2 + b^2] - a]/Sqrt[2] for a,b > 0

or

0 < b < 2*c*Sqrt[a + c^2] for a,c > 0

or

a > (b^2 - 4*c^4)/(4*c^2) for b,c > 0

are met. I hope I got the conditions right...

Edit:

actually one condition is enough:

b < 2*c*Sqrt[a + c^2]

With it the Log replacements are correct and you get:

Simplify[Simplify[Together[
ToRadicals[
Integrate[1/(a^2 + b^2 + 4*b*c*x + 2*a*(c^2 - x^2) + 
                (c^2 + x^2)^2), {x, -Infinity, 0}] + 
Integrate[1/(a^2 + b^2 + 4*b*c*x + 2*a*(c^2 - x^2) + 
(c^2 + x^2)^2),{x, 0, Infinity}]]]] /. 
    {Log[-Sqrt[a - I*b] - I*c] -> 
     Log[Sqrt[a - I*b] + I*c] - I*Pi, 
     Log[-Sqrt[a + I*b] + I*c] -> 
     Log[Sqrt[a + I*b] - I*c] + I*Pi, 
     Log[-Sqrt[a - I*b] + I*c] -> 
     Log[Sqrt[a - I*b] - I*c] + I*Pi, 
     Log[-Sqrt[a + I*b] - I*c] -> 
     Log[Sqrt[a + I*b] + I*c] - I*Pi}]
(* (2*c*Pi)/(-b^2 + 4*c^2*(a + c^2)) *)

Tested with plots (following code was edited + condition for range of c added):

a = RandomReal[{0, 3}]; 
b = RandomReal[{-3, 3}]; 
c = RandomReal[{Sqrt[Abs[b]/2], 3}]; 
p1 = LogPlot[{NIntegrate[1/(a^2 + b^2 + 4*b*c*x + 
 2*a*(c^2 - x^2) + (c^2 + x^2)^2), {x, -100, 100}, 
 WorkingPrecision -> 200], (2*Pi*c)/(4*c^2*(a + c^2) - b^2)}, 
 {a, (b^2 - 4*c^4)/(4*c^2), 8}, PlotStyle -> {Blue, Dashed}]; 
p2 = LogPlot[{NIntegrate[1/(a^2 + b^2 + 4*b*c*x + 
 2*a*(c^2 - x^2) + (c^2 + x^2)^2), {x, -100, 100}, 
 WorkingPrecision -> 200], (2*Pi*c)/(4*c^2*(a + c^2) - b^2)}, 
 {b, -2*c*Sqrt[a + c^2], 2*c*Sqrt[a + c^2]}, PlotStyle -> {Green, Dashed}]; 
p3 = LogPlot[{NIntegrate[1/(a^2 + b^2 + 4*b*c*x + 
 2*a*(c^2 - x^2) + (c^2 + x^2)^2), {x, -100, 100}, 
 WorkingPrecision -> 200], (2*Pi*c)/(4*c^2*(a + c^2) - b^2)}, 
 {c, Sqrt[Sqrt[a^2 + b^2] - a]/Sqrt[2], 8}, PlotStyle -> {Red, Dashed}]; 
{a, b, c}
Show[p1, p2, p3, PlotRange -> All]

enter image description here

Improve this answer
$\endgroup$
13
  • $\begingroup$ Can you elaborate " Split the integration range from -Infinity to 0 and from 0 to Infinity. Add the two Rootsums and convert the result to radicals, then help Mathematica simplify the result"? In other case, all that is built on the sand. $\endgroup$
    – user64494
    Commented Mar 26, 2022 at 15:38
  • $\begingroup$ @user64494 Well, over the whole range the integral reports zero for some reason. the split gives two rootsums which then have to be added, apply ToRadicals[] and simplify $\endgroup$
    – Andreas
    Commented Mar 26, 2022 at 17:33
  • 1
    $\begingroup$ Sorry, I am waiting a Mathematica code, not empty words. $\endgroup$
    – user64494
    Commented Mar 26, 2022 at 17:42
  • $\begingroup$ @user64494 don't pretend you have never used Mathematica before. I wrote everything necessary to test the calculation and won't expand on this more $\endgroup$
    – Andreas
    Commented Mar 26, 2022 at 17:47
  • 3
    $\begingroup$ @Andreas I would add that, for me, trying Integrate[f, {x, -Infinity, 0, Infinity},...] is a standard move when Integrate misbehaves. I'm not sure why putting a waypoint on the integration path sometimes helps, but I think it changes how Integrate deals with branch points. $\endgroup$
    – Michael E2
    Commented Mar 30, 2022 at 15:17
5
$\begingroup$

Edit2 corrected error in cc2 and int2

Now for some paramter values int1 is valid, for other int2 and for some int1 yield wrong complex solutions. See Manipulate...

Edit Of course straightforward integration should give the right result. But it does not. This is a bug !

Nevertheless found a way to elict the right result, but still some partly wrong answer in conditions remain.

Split the integrand with Apart.

f[a_, b_, c_] = 
    1/((-(a - I b) + (x + I c)^2) (-(a + I b) + (x - I c)^2));

fap = f[a, b, c] // Apart

(*   (-I b + 4 c^2 + 2 I c x)/(
 2 (b^2 - 4 a c^2 - 4 c^4) (-a - I b - c^2 - 2 I c x + x^2)) + (
 I b + 4 c^2 - 2 I c x)/(
 2 (b^2 - 4 a c^2 - 4 c^4) (-a + I b - c^2 + 2 I c x + x^2))   *)

cc1 = ComplexExpand[fap[[1]], TargetFunctions -> {Re, Im}] // 
  Simplify[#, Assumptions -> c > 0 && a > 0 && Element[b, Reals]] &

(*   (I (b + 4 I c^2 - 2 c x))/(2 (b^2 - 4 c^2 (a + c^2)) (a + 
   I b + (c + I x)^2))   *)

cc2 = ComplexExpand[fap[[2]], TargetFunctions -> {Re, Im}] // 
  Simplify[#, Assumptions -> c > 0 && a > 0 && Element[b, Reals]] &

    (*   (b - 4 I c^2 - 
 2 c x)/(2 (b^2 - 4 c^2 (a + c^2)) (I a + b + I (c - I x)^2))   *)


int1[a_, b_, c_] = 
 Integrate[cc1, {x, -\[Infinity], \[Infinity]}, 
   Assumptions -> c > 0 && a > 0 && Element[b, Reals]] // 
  FullSimplify[#, 
    Assumptions -> c > 0 && a > 0 && b \[Element] Reals] &

(*   ConditionalExpression[-((2 c \[Pi])/(
  b^2 - 4 c^2 (a + c^2))), -c - Im[Sqrt[a + I b]] < 
   0 && -c + Im[Sqrt[a + I b]] < 
   0 && (b <= 2 c^2 || 4 c^2 (a + c^2) > b^2)]   *)
int2[a_, b_, c_] = 
 Integrate[cc2, {x, -\[Infinity], \[Infinity]}, 
   Assumptions -> c > 0 && a > 0 && Element[b, Reals]] // 
  FullSimplify[#, 
    Assumptions -> c > 0 && a > 0 && b \[Element] Reals] &

(*   ConditionalExpression[((b + 2 (Sqrt[a - I b] - I c) c) \[Pi])/(
 2 Sqrt[a - I b] (b^2 - 4 c^2 (a + c^2))), 
 b > 2 c^2 && 4 c^2 (a + c^2) < b^2]   *)

int1[1, 1, 1] // N
(*   0.897598   *)

nint[a_, b_, c_] := 
   NIntegrate[f[a, b, c], {x, -\[Infinity], \[Infinity]}]

nint[1, 1, 1]
(*   0.897598\[VeryThinSpace]+ 0. I   *)

Manipulate[{int1[a, b, c] // N, int2[a, b, c] // N, 
  nint[a, b, c] // Chop}, {{a, 1}, 0, 5}, {{b, 1}, -5, 5}, {{c, 1}, 0,
   6}]

$Version  (*   "8.0 for Microsoft Windows (32-bit) (December 9, 2010)"   *)

Edit3

Show, where int1 and int2 are defined.

RegionPlot3D[-c - Im[Sqrt[a + I b]] < 0 && -c + Im[Sqrt[a + I b]] < 
   0 && (b <= 2 c^2 || 4 c^2 (a + c^2) > b^2), {a, 0, 10}, {b, -5, 
  5}, {c, 0, 10}]

RegionPlot3D[
 b > 2 c^2 && 4 c^2 (a + c^2) < b^2, {a, 0, 10}, {b, -5, 5}, {c, 0, 
  10}]
Improve this answer
$\endgroup$
8
  • $\begingroup$ -1. Unfortunately, Im[Root[a^2 + b^2 + 2 a c^2 + c^4 + 4 b c #1 + (-2 a + 2 c^2) #1^2 + #1^4 &, 3]] > 0 && Im[Root[a^2 + b^2 + 2 a c^2 + c^4 + 4 b c #1 + (-2 a + 2 c^2) #1^2 + #1^4 &, 2]] > 0 && Im[Root[a^2 + b^2 + 2 a c^2 + c^4 + 4 b c #1 + (-2 a + 2 c^2) #1^2 + #1^4 &, 4]] > 0 && Im[Root[a^2 + b^2 + 2 a c^2 + c^4 + 4 b c #1 + (-2 a + 2 c^2) #1^2 + #1^4 &, 1]] > 0 /. {a -> 1, b -> 1, c -> 1} produces False. Therefore, int2[1, 1, 1] // N is wrong. $\endgroup$
    – user64494
    Commented Mar 26, 2022 at 15:35
  • $\begingroup$ Also cannot execute f[a_, b_, c_] = 1/((-(a - I b) + (x + I c)^2) (-(a + I b) + (x - I c)^2)); fap = f[a, b, c] // Apart; cc1 = ComplexExpand[fap[[1]], TargetFunctions -> {Re, Im}] // Simplify[#, Assumptions -> c > 0 && a > 0 && Element[b, Reals]] &, obtaining "Tag Times in (1/((-a-I\b+(-I c+x)^2)(-a+I\b+(I c+x)^2)))[a_,b_,c_] \ is Protected". $\endgroup$
    – user64494
    Commented Mar 26, 2022 at 15:47
  • $\begingroup$ ClearAll helps for the above. But executing int1[a_, b_, c_] = Integrate[cc1, {x, -\[Infinity], \[Infinity]}, Assumptions -> c > 0 && a > 0 && Element[b, Reals]] // FullSimplify[#, Assumptions -> c > 0 && a > 0 && b \[Element] Reals] &, I obtain "Integral of (I (b+4 I c^2-2 c x))/(2 (b^2-4 c^2 (a+c^2)) (a+I b+(c+I x)^2)) does not converge on {-[Infinity],[Infinity]" and the returned input in version 13 on Windows 10. $\endgroup$
    – user64494
    Commented Mar 26, 2022 at 17:38
  • $\begingroup$ Please present the executed code as nb. file and/or its pdf file through Dropbox or other tool. In opposite case your answer is empty words (Version 8 is rare at the present.). Hope I am clear. Regard. $\endgroup$
    – user64494
    Commented Mar 26, 2022 at 18:18
  • $\begingroup$ You wrote "This is a bug !". This is not a bug, but a weakness. Version 13 on Windows produces 0 under an impossible condition. $\endgroup$
    – user64494
    Commented Mar 26, 2022 at 18:23
-6
$\begingroup$

Mathematica is right. Indeed, the command

Integrate[ 1/((-(a - I b) + (x + I c)^2) (-(a + I b) + (x - I c)^2)),
 {x, -Infinity, Infinity}, Assumptions -> a > 0 && b > 0 && c \[Element] Reals]

ConditionalExpression[0, Im[Root[a^2 + b^2 + 2*a*c^2 + c^4 + 4*b*c*#1 + (-2*a + 2*c^2)*#1^2 + #1^4 & , 1]] > 0 && Im[Root[a^2 + b^2 + 2*a*c^2 + c^4 + 4*b*c*#1 + (-2*a + 2*c^2)*#1^2 + #1^4 & , 2]] > 0 && Im[Root[a^2 + b^2 + 2*a*c^2 + c^4 + 4*b*c*#1 + (-2*a + 2*c^2)*#1^2 + #1^4 & , 3]] > 0 && Im[Root[a^2 + b^2 + 2*a*c^2 + c^4 + 4*b*c*#1 + (-2*a + 2*c^2)*#1^2 + #1^4 & , 4]] > 0]

results in 0, but the condition is impossible: all the four roots of a polynomial of fourth degree with real coefficients cannot have the positive imaginary parts.

PS. So Mathematica fails here, not producing a constructive answer. However, its answer is not wrong: under the certain impossible condition the integral equals 0.

Edit. Skip of key-board: "imaginary" is added.

Improve this answer
$\endgroup$
7
  • $\begingroup$ Please explain the down-vote. What is incorrect in my answer? $\endgroup$
    – user64494
    Commented Mar 26, 2022 at 12:36
  • $\begingroup$ –1: your answer is obviously wrong. A simple With[{a = 0.2, b = 0.3, c = 0.4}, LogPlot[1/((-(a - I b) + (x + I c)^2) (-(a + I b) + (x - I c)^2)), {x, -10, 10}]] shows that the integrand is positive-definite and therefore the integral must be positive for these specific values of $a, b, c$ (which I've picked at random). $\endgroup$
    – Roman
    Commented Mar 28, 2022 at 8:16
  • $\begingroup$ @Roman: Master your math. The Mathematica answer claims that the integral equals zero under the certain impossible condition. Don't hesitate to ask for further explanation in need. $\endgroup$
    – user64494
    Commented Mar 30, 2022 at 5:23
  • 3
    $\begingroup$ Technically you are right; in practice, this answer is useless because wrong for all values of $a,b,c$. What the OP wants is an expression that gives the integral in question, and not a counter-factual but "technically correct" mathematical subtlety. "The integral is zero whenever your parameters satisfy a non-satisfiable condition" is not a useful answer; calling it "mathematically correct" does not help in practice. $\endgroup$
    – Roman
    Commented Mar 30, 2022 at 8:00
  • 1
    $\begingroup$ @Roman has hit the heart of the matter: mathematics as practiced by most of us is useful. Sophistry is useless. $\endgroup$
    – John Doty
    Commented Mar 30, 2022 at 12:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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