# Zero result from a positive integral

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.

• 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. Mar 26, 2022 at 11:28
• @user64494 The integrand is positive.... Mar 26, 2022 at 11:43
• xiaohuamao (@ does not work): The integrand is complex-valued, not positive. Master your math. Mar 26, 2022 at 12:10
• @user64494 The integrand is in the form $a\times a^*$, which is positive. Mar 26, 2022 at 12:16
• xiaohuamao (@ does not work): You are right, Mathematica too. Mar 26, 2022 at 12:32

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[
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]


• 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. Mar 26, 2022 at 15:38
• @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 Mar 26, 2022 at 17:33
• Sorry, I am waiting a Mathematica code, not empty words. Mar 26, 2022 at 17:42
• @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 Mar 26, 2022 at 17:47
• @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. Mar 30, 2022 at 15:17

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}]  • -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. Mar 26, 2022 at 15:35 • 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". Mar 26, 2022 at 15:47 • 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. Mar 26, 2022 at 17:38 • 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. Mar 26, 2022 at 18:18 • You wrote "This is a bug !". This is not a bug, but a weakness. Version 13 on Windows produces 0 under an impossible condition. Mar 26, 2022 at 18:23 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. • Please explain the down-vote. What is incorrect in my answer? Mar 26, 2022 at 12:36 • –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). Mar 28, 2022 at 8:16 • @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. Mar 30, 2022 at 5:23 • 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. Mar 30, 2022 at 8:00
• @Roman has hit the heart of the matter: mathematics as practiced by most of us is useful. Sophistry is useless. Mar 30, 2022 at 12:44