Why this real integral yields imaginary results?

This integral yields -1-4Iπ/3 in Mathematica:

Integrate[(y - y^2 + x - x^2 + 2*x*y)/(1 - x - y), {x,0,1}, {y, 0, 1}]

Since the integrand is real (although divergent at some points), how does Mathematica come up with a complex number?

• I think it's a bug. Have you reported this to WRI? – xzczd Sep 27 '16 at 5:49
• I don't think it's a bug, since sympy and other programs give the same result. – Chong Wang Sep 27 '16 at 5:55
• I think it's just an possible evidence that Sympy has used the same algorithm as Mathematica :). See my answer below. – xzczd Sep 27 '16 at 6:09
• After a second thought, I noticed it's indeed a bug. The integration is divergent. See my update. – xzczd Sep 27 '16 at 9:17
• (at) Chong Wang: The integral is obviously divergent. Hence it must be DEFINED first of all. bbgodfrey has given n possibility in his answer. – Dr. Wolfgang Hintze Sep 27 '16 at 9:27

To assure that the path of integration stays on the real axis in x and y, use

Integrate[(y - y^2 + x - x^2 + 2*x*y)/(1 - x - y), {x, 0, 1}, {y, 0, 1},
PrincipalValue -> True]

(* -1 *)

Alternative approach

Another approach, more complicated but perhaps informative, is to solve the integral using coordinates, {u == x + y, v == x - y}.

Flatten@Solve[{u == x + y, v == x - y}, {x, y}]
Simplify[(y - y^2 + x - x^2 + 2*x*y)/(1 - x - y)/2 /. tr]
(* (u - v^2)/(2 - 2 u) *)

(Division by two above is needed to take account of the Jacobian of the transformation.) The integral over v does not involve the singularity and can be performed without difficulty for both u < 1 and u > 1.

um = Simplify@Integrate[%, {v, -u, u}] // Apart
(* -(2/3) - 2/(3 (-1 + u)) - (2 u)/3 + u^2/3 *)
up = Simplify@Integrate[%%, {v, -2 + u, 2 - u}] // Apart
(* -(10/3) - 2/(3 (-1 + u)) + (8 u)/3 - u^2/3 *)

The singular term is identical for both expressions and can be separated to yield

sv1 = Piecewise[{{um + 2/(3 (-1 + u)), u <= 1}, {up + 2/(3 (-1 + u)), u > 1}}];
Integrate[sv1, {u, 0, 2}] -
Integrate[2/(3 (-1 + u)), {u, 0, 2}, PrincipalValue -> True]
(* -1 *)

(The second integral is identically zero.) For completeness,

Plot[sv1, {u, 0, 2}] Note that Integrate has difficulties, if sv1 and 2/(3 (-1 + u)) are included in the same integrand, and I posed question 127448 to seek advice on why. xzczd provided a good workaround.

• What will Mathematica do if I don't insist principle value? – Chong Wang Sep 27 '16 at 5:49
• @ChongWang It probably does a complex contour integral to avoid the singularity. At least that would explain the complex result. – Michael E2 Sep 27 '16 at 10:36
• Typo?: I think you dropped a minus sign. I get -1. – Michael E2 Sep 27 '16 at 11:16
• Hopefully anyone using this will be sure they understand what principal value means. (i.e. read Dr Hintz answer). – george2079 Sep 27 '16 at 17:55
• @Myridium I think Mathematica treats the integral as an iterated integral, so each 1D integral, the interior one over y depending on x as a parameter, might be approached as a line integral in the complex plane. By the Cauchy integral formula, you cannot get any result you like, but only ones that differ by multiples of -4 I Pi/3, depending how the interior integral winds around the pole at y == 1 - x. (I'm assuming a path that does not depend on x and does not pass through a pole for any x.) – Michael E2 Sep 28 '16 at 10:22

Recently, in this forum divergent integrals have become kind of abundant. Therefore, it might be worth to dwell a little more on the subject.

Abstract

The integral of the OP is obviously divergent. Hence, in order to give it sense at all, the divergence must be circumvented by some additional prescription, normally called "regularization".

We study two types of regularization procedures, taking a principal value and using a parametric integral.

We shall see that the final result depends on the regularization procedure. So that, strictly speaking, the integral can have any value we want.

Not issuing an error message stating the divergence of the integral is considered a deficiency of Mathematica.

Restatement of the problem

With the integrand

h = (y - y^2 + x - x^2 + 2*x*y)/(1 - x - y);

the integral in question is

f0 := Integrate[h, {x, 0, 1}, {y, 0, 1}]

The integral is calculated by Mathematica in a few seconds, giving.

f0

(* Out= -1 - (4 I π)/3 *)

Is this really wrong? Can an integral over a real integrand give a complex result? Yes, it can, because the integral is divergent.

Notice that Mathematica does not give the error message that the integral in divergent, what it certainly is.

This missing error message a deficiency which could be called a bug.

Let's now turn to the regularization.

Principal value (p.v.)

The problem of (certan types of) divergent integrals has been studied already by the famous Cauchy who recommended to take his principal value. This has been done here in the answer of bbgoodfrey.

Now Cauchy's principal value implemented in Mathematica as an option PrincipalValue->True to Integrate[] is a special case of a symmteric approach from both sides of the singularity.

There are others possibilities, however.

Consider the (divergent) integral

f = Integrate[1/x, {x, -1, 1}]

> During evaluation of In:= Integrate::idiv: Integral of 1/x does
> not converge on {-1,1}. >>

(* Out= \!$$\*SubsuperscriptBox[\(\[Integral]$$, $$-1$$, $$1$$]$$\*FractionBox[\(1$$, $$x$$] \[DifferentialD]x\)\) *)

Here Mathematica correctly complains that the integral is divergent.

Cauchy's p.v. gives

fc = Integrate[1/x, {x, -1, 1}, PrincipalValue -> True]

(* Out= 0 *)

Now take the general p.v. which is not necessarily symmetric but defined by

Limit[Integrate[1/x, {x, -1, -a ϵ},
Assumptions -> {a > 0, b > 0, ϵ > 0}] +
Integrate[1/x, {x, b ϵ, 1},
Assumptions -> {a > 0, b > 0, ϵ > 0}], ϵ -> 0]

(* Out= Log[a] - Log[b] *)

Hence the integral can attain any real value depending of the choice of a and b. Only if a == b the p.v. is zero (Cauchy case).

Now let us study the general p.v. in the example of the two-dimensional integral of the OP.

This can be done using Boole[].

For the lower triangle in the xy-region the integral is

lower = Integrate[(y - y^2 + x - x^2 + 2*x*y)/(1 - x - y)
Boole[x + y <= 1 - a ϵ], {x, 0, 1}, {y, 0, 1},
Assumptions -> {1/10 > a > 0, 0 < ϵ < 1/100}]

(* Out= 1/9 (-8 + 9 a ϵ - a^3 ϵ^3 - 6 Log[a ϵ]) *)

whereas in the upper triangle we have

upper = Integrate[(y - y^2 + x - x^2 + 2*x*y)/(1 - x - y)
Boole[x + y >= 1 + b ϵ], {x, 0, 1}, {y, 0, 1},
Assumptions -> {1/10 > b > 0, 0 < ϵ < 1/100}]

(* Out= 1/9 (-1 + 9 b ϵ - 9 b^2 ϵ^2 + b^3 ϵ^3 +
6 Log[b ϵ]) *)

The sum is

all = lower + upper // Simplify

(* Out= 1/9 (-9 + 9 a ϵ + 9 b ϵ - 9 b^2 ϵ^2 -
a^3 ϵ^3 + b^3 ϵ^3 - 6 Log[a ϵ] +
6 Log[b ϵ]) *)

And the regularized integral becomes the limit ϵ->0 of this

fpv = Limit[all, ϵ -> 0] // Expand

(* Out= -1 - (2 Log[a])/3 + (2 Log[b])/3 *)

Hence, beside the "well known" -1 we get an arbitrary additional real summand.

Parametric integral

We can get rid of the "annoying denominator" using the representation

Integrate[Exp[ t (1 - x - y)], t] /. t -> 0

(* Out= -(1/(-1 + x + y)) *)

Our integral then becomes dependent on the parameter t. Mathematica calculates the x-integral explicitly

fpx = Integrate[(y - y^2 + x - x^2 + 2*x*y) Exp[ t (1 - x - y)], {x, 0,
1}, {y, 0, 1}, Assumptions -> t > 0]

(* Out= (2 (2 - 2 t + (-2 + t (2 + t)) Cosh[t] - t^2 Sinh[t]))/t^4 *)

Now the indefinite t-integral becomes

fpt = Integrate[fpx, t]

(* Out= (-4 + 6 t + (4 - 6 t - 4 t^2) Cosh[t] + 2 t Sinh[t] +
4 t^3 SinhIntegral[t])/(3 t^3) *)

Taking the now the limit we get

Limit[fpt, t -> 0]

(* Out= -1 *)

Summary

Two of the considered regularization methos yield -1 as the result. Only the somewhat "exotic" unsymmteric p.v. gives an arbitrary real result.

EDIT #1

The critical remark towards Mathamatica can be softened by observing that doing only one integral (x) of the two and providing the range of the other variable (y) as an assumption we get

Integrate[h, {x, 0, 1}, Assumptions -> 0 < y < 1]

During evaluation of In:= Integrate::idiv: Integral of (x+y-(-x+y)^2)/(-1+x+y) does not converge on {0,1}. >>

(* Out= Integrate[(x - (x - y)^2 + y)/(1 - x - y), {x, 0, 1},
Assumptions -> 0 < y < 1] *)

Here we have the correct error message and the unevaluated integral.

EDIT #2

We can also regularize by analytic continuation. This will lead to complex results.

Replacing the 1 in the denominator by a variable t and calculating the integral in a region of t where the integral is convergent gives

ft = Integrate[(x + y - (x - y)^2)/(t - x - y), {x, 0, 1}, {y, 0, 1},
Assumptions -> t > 2]

(* Out= 1/3 (-8 + 5 t +
4 (4 + 3 (-3 + t) t) ArcCoth[
3 - 2 t] - (-3 + t) t^2 (Log[-2 + t] - 2 Log[-1 + t] + Log[t])) *)

This function has a branch cut fromt = 0 to t = 2. This includes the point t = 1. Hence we can expext complex values as well as more than one value.

On the real t-axis we have graphically

Plot[{Re[ft], Im[ft]}, {t, -2, 5}, PlotRange -> All,
PlotLabel ->
"The function ft(t)\nRe[ft] - blue curve\nIm[ft] - yellow curve",
AxesLabel -> {"t", "Re,Im(f)"},
Epilog -> {Red, PointSize -> 0.02, Point[{1, 0}]}] The limit of ft for t->1 from below is

Limit[ft, t -> 1, Direction -> +1]

(* Out= -1 - (2 I \[Pi])/3 *)

Ths values is the same as with the other regularizations discussed.

From above we find an new value

Limit[ft, t -> 1, Direction -> -1]

(* Out= -1 + 2 I \[Pi] *)

Hence the real part is -1 as expected, and the jumps show that we have two different imaginary parts at t = 1.

The complete picture of ft for complex t is

Plot3D[Re[ft] /. t -> u + I v, {u, 0, 2}, {v, -1/2, 1/2},
PlotLabel -> "Re[ft]", AxesLabel -> {"Re[t]", "Im[t]", "Re[ft]"}] Plot3D[Im[ft] /. t -> u + I v, {u, -1/2, 2 + 1/2}, {v, -1, 1},
PlotLabel -> "Im[ft]", AxesLabel -> {"Re[t]", "Im[t]", "Im[ft]"},
PlotRange -> All] Unfortunately, I don't know how to rotate the output of Plot3D so that the cut would be better visible. I'm sure someone else can help me.

PS: I am curious in which manner Michael E4 was able to produce even more different imaginary parts from the integral.

EDIT #3 Simplifcation of the problem

The question was: "how does Mathematica come up with a complex number?"

Let us simplify the problem down to the core and consider this integral

g = Integrate[1/(1 - x - y), {x, 0, 1}, {y, 0, 1}]

Out= I \[Pi]

We do the integral step by step.

First the x-integral:

Integrate[1/(1 - x - y), {x, 0, 1}]

Out= ConditionalExpression[Log[-1 + y] - Log[y],
Re[y] > 1 || Re[y] < 0 || y \[NotElement] Reals]

We select a condition and repeat the integration

Integrate[1/(1 - x - y), {x, 0, 1}, Assumptions -> y \[NotElement] Reals]

Out= Log[-1 + y] - Log[y]

Since the argument of the first Log is negative in the range of interest we make use of the relation

Simplify[Log[-1 + y] == I \[Pi] + Log[1 - y], 0 < y < 1]

Out= True

This transformation leads to a complex number.

The y-integral is then finally

Integrate[I \[Pi] + Log[1 - y] - Log[y], {y, 0, 1}]

Out= I \[Pi]

The complex number result was produced making contradictory assumptions with respect to y being real or not real. Hence the result can be considered as spurious. On the other hand remember that the integral is divergent. So we can consider the development here boldly as a regularization.

The principal value of the integral is zero

gp = Integrate[1/(1 - x - y), {x, 0, 1}, {y, 0, 1}, PrincipalValue -> True]

Out= 0

Another simple integral gives a complex number with non-zero real part

g2 = Integrate[x/(1 - x - y), {x, 0, 1}, {y, 0, 1}]

Out= 1/2 I (I + \[Pi])

Expand[%]

Out= -(1/2) + (I \[Pi])/2

The principal values gives the real part of the integral

g2p = Integrate[x/(1 - x - y), {x, 0, 1}, {y, 0, 1}, PrincipalValue -> True]

Out= -(1/2)
• Re edit #1: This is effectively what Mathematica does with the interior integral when one sets GenerateConditions -> True, since True is the default for single integrals and it automatically assumes the integration range for the outer variable. – Michael E2 Sep 27 '16 at 12:22
• I'm giving up for now on figuring out the -4 I π/3. I can get -1 + 2 I π/3 + 4 I π n/3 for any integer n. The imaginary part may be a bug. Any ideas? – Michael E2 Sep 27 '16 at 13:47
• @Michael E4: I can get only -4 I π/3 and -1 - 2 I π/3 – Dr. Wolfgang Hintze Sep 27 '16 at 18:40
• @Michael E4 see my EDIT #3 for a simplified version of the problem. – Dr. Wolfgang Hintze Sep 29 '16 at 7:32

In multiple integrals, GenerateConditions -> Automatic is set to False in interior integrals (What exactly does GenerateConditions do?).

If we override it, Mathematica reports that the integral is divergent.

Integrate[(y - y^2 + x - x^2 + 2*x*y)/(1 - x - y),
{x, 0, 1}, {y, 0, 1}, GenerateConditions -> True]

Integrate::idiv: Integral of -(x/(-1+x+y))+x^2/(-1+x+y)-y/(-1+x+y)-(2 x y)/(-1+x+y)+y^2/(-1+x+y) does not converge on {0,1}.

I would say it's not a bug, since this behavior seems intended. There is a balancing act in performance between speed and accuracy that GenerateConditions -> Automatic addresses. However, this behavior should probably be made more widely known (say, by putting it in the docs for Integrate). A complex result for a real integral may indicate something went awry, but it's not as helpful a hint as warning such as "Complex result for a real integral: Integral may be divergent, etc.; try GenerateConditions -> True."

Whether PrincipalValue -> True is a workaround depends on the application. Sometimes a divergent integral is simply divergent.

• Interestingly, NIntegrate returns close to -1, suggesting that the integral does converge; but doing NIntegrate on both halves of the integral ($0<y<1-x$ and $1-x<y<1$) shows that each half (and the whole integral) is actually divergent. – 2012rcampion Sep 27 '16 at 16:25
• @2012rcampion The approximate -1 is due in part to the symmetry of the integrand (int - (int /. {x -> 1 - x, y -> 1 - y}) // Simplify), so that the values at symmetric sample points average to -1. (The sampling seems to be imperfectly symmetric.) As you ramp up MaxRecursion the result becomes unstable, due rounding error from large integrand values near the singularity. (Increasing WorkingPrecision can then control the rounding error, too.) A little asymmetry, such as {y, 0, 1 - $MachineEpsilon}, throws the value way off, like Dr. Hintze's asymmetric p.v. – Michael E2 Sep 27 '16 at 16:55 • Yes (and upvote), this default behavior of convergence checking is by intent, and for the reasons given above. – Daniel Lichtblau Sep 28 '16 at 0:22 Update After a second thought, I noticed this is indeed a bug, and the integration is divergent. Just recall the truth that Integrate[1/x, {x, -1, 1}] is divergent, and notice OP's integration in essential belongs to the same class. Actually if we rotate the integrand for 45 degree, Mathematica will successfully notice the integration is divergent: rotatedintegrand = (y - y^2 + x - x^2 + 2 x y)/(1 - x - y) /. {x -> (ξ + η)/Sqrt, y -> (ξ - η)/Sqrt} // Simplify (* (2 η^2 - Sqrt ξ)/(-1 + Sqrt ξ) *) excludedarea = Abs[ξ - Sqrt/2] < Sqrt/2 && Abs@η < Sqrt/2 && Abs[ξ - Sqrt/2] + Abs[η] > Sqrt/2; newint = Integrate[#, {ξ, 0, Sqrt}, {η, -Sqrt/2, Sqrt/2}] &; newint[rotatedintegrand] - newint[Boole@excludedarea rotatedintegrand] Integrate::idiv Original Answer As mentioned in the comment above, I think it's a bug, and one work-around I found is: Integrate[(y - y^2 + x - x^2 + 2 x y)/(1 - x - y) /. x -> x + 2, {x, -2, -1}, {y, 0, 1}] (* -1 *) Notice I've only done a trivial shift on the integrand, the integration is equivalent to the original. The offset -2 is found by trial and error, using other offset may lead to other incorrect result, for example: Integrate[(y - y^2 + x - x^2 + 2 x y)/(1 - x - y) /. x -> x + 3, {x, -3, -2}, {y, 0, 1}] (* -1 + (I π)/6 *) • This is really weired. I guess I need more time to read the documentation. But how to report it to WRI? – Chong Wang Sep 27 '16 at 6:51 • @ChongWang You can report possible bug on this page: wolfram.com/support/contact/email – xzczd Sep 27 '16 at 6:53 • BTW, to the downvoter, I am interested in what was missing from my answer. Did I give too little detail? Or, was it considered incorrect in some way? Either way, would you please elaborate? I'm not trying to complain here, I'm just curious about what I could have done instead. – xzczd Sep 27 '16 at 6:55 • -1: Considering some behaviour a Bug is not an answer. Your workaround may solve the case, but gives no answer to the Headline-Question. Additionally, @bbgodfrey found an "explanation-like" answer, which led me to the interesting read of ref/PrincipalValue: "Setting PrincipalValue->True gives finite answers for integrals that had simple pole divergences with PrincipalValue->False." – DPF Sep 27 '16 at 7:09 • Your additional explanation of your thoughts lets me draw back the downvote. – DPF Sep 27 '16 at 7:41 In case this might shed more light, I used Maxima, which asks the user questions to decide which branch to take. This is a summary of what the choices it asked about. All are related to the range of$y$and not$x\$. In case I made an error in the table above, here is the actual maxima session Noticed also that Maxima gives -1 - (2 I Pi)/3 while Mathematica gives -1 - (4 I Pi)/3 for some of these values.

Maple gives undefined on this integral. When trying to use Assumptions in Mathematica to to try to duplicate Maxima result, it did not give same values. Here is my attempt:

Integrate[(y - y^2 + x - x^2 + 2*x*y)/(1 - x - y), {x, 0, 1}, {y, 0,1},
Assumptions -> {y > 1}]

etc... Summary of the above: 