# Real-valued integral gives large imaginary part

I found a very strange result in Integrate command. I numerically integrated $$\int_{-1}^1 \int_{-1}^1 \frac{1}{|x|^{2/5} |y|^{3/10} |x-y|^{7/10}} dx dy.$$

int = 1/(Abs[x]^(2/5) Abs[y]^(3/10) Abs[x - y]^(7/10))
val = Integrate[int, {x, -1, 1}, {y, -1, 1}]
N[val]


gives the complex answer 17.0597 - 13.2162 I ! (Before evaluating numerically, the Mathematica gives an analytic expression for val as

(147 2^(4/5) Sqrt[\[Pi]] Gamma[-(3/10)] Gamma[1/5] -
210 \[Pi] Gamma[-(3/10)] Gamma[17/10] +
210 Sqrt[5] \[Pi] Gamma[-(3/10)] Gamma[17/10] -
700 (-2)^(4/5) Sqrt[\[Pi]] Gamma[1/5] Gamma[17/10] -
2100 Gamma[-(3/10)] Gamma[17/10] Hypergeometric2F1[1/10, 7/10, 11/
10, -1] +
300 Gamma[-(3/10)] Gamma[17/10] Hypergeometric2F1[7/10, 7/10, 17/
10, -1] -
700 Sqrt[5 + Sqrt[5]] Gamma[-(3/5)] Gamma[13/10] Gamma[17/10] Root[
512 + #^10& , 4, 0])/(63 Gamma[-(3/10)] Gamma[17/10])


.) Note that the large imaginary part cannot be regarded as the numerical error. This is very strange, since the integral is real-valued. Why this error occurs?

• I am using version 12.1.
• The same issue in 12.3.1 on Windows 10. RealAbs instead of Abs does not help. Oct 13 at 9:35

This seems to be a bug in higher versions of MMA.

In version "8.0 for Microsoft Windows (32-bit) (December 9, 2010)"  you get the right result independent of the order of integration. Although the two soutions look different, they are equivalent.

integrand = 1/(Abs[x]^(2/5) Abs[y]^(3/10) Abs[x - y]^(7/10));

int1 = Integrate[integrand, {x, -1, 1}, {y, -1, 1}]

int2 = Integrate[integrand, {y, -1, 1}, {x, -1, 1}]

(*   (49 \[Pi]^(3/2) Gamma[1/10] Gamma[7/10] Gamma[9/10] -
49 Sqrt[5] \[Pi]^(3/2) Gamma[1/10] Gamma[7/10] Gamma[9/10] -
49 \[Pi]^(3/2) Gamma[2/5] Gamma[7/10]^2 Gamma[9/10] +
49 Sqrt[5] \[Pi]^(3/2) Gamma[2/5] Gamma[7/10]^2 Gamma[9/10] +
70 Sqrt[\[Pi]]
Gamma[3/10] Gamma[2/5] Gamma[3/5] Gamma[7/10] Gamma[17/10] +
35 2^(4/5) \[Pi] Gamma[1/5] Gamma[2/5] Gamma[9/10] Gamma[17/10] +
35 2^(3/10) Sqrt[5 + Sqrt[5]]
Gamma[1/5] Gamma[3/10] Gamma[2/5] Gamma[7/10] Gamma[9/10] Gamma[
17/10] +
175 2^(4/5) \[Pi] Gamma[2/5] Gamma[9/10] Gamma[6/5] Gamma[17/10] -
700 Sqrt[\[Pi]]
Gamma[2/5] Gamma[7/10] Gamma[9/10] Gamma[17/
10] Hypergeometric2F1[1/10, 7/10, 11/10, -1] +
100 Sqrt[\[Pi]]
Gamma[2/5] Gamma[7/10] Gamma[9/10] Gamma[17/
10] Hypergeometric2F1[7/10, 7/10, 17/10, -1])/(21 Sqrt[\[Pi]]
Gamma[2/5] Gamma[7/10] Gamma[9/10] Gamma[17/10])   *)

int2 // N

(*   35.2503   *)

int1 - int2 // FullSimplify

(*   0   *)


Edit

May be other versions can do the job, if you get rid of the Abs. Split integration region into 3 parts according to plot of integrand. (Factor 2 due to symmetry)

Plot3D[integrand, {x, -1, 1}, {y, -1, 1}, PlotPoints -> 60]

i1 = PowerExpand[integrand, Assumptions -> {0 < x < 1, -1 < y < 0}]

i2 = PowerExpand[integrand, Assumptions -> {0 < x < 1, 0 < y < x}]

i3 = PowerExpand[integrand, Assumptions -> {0 < x < 1, x < y < 1}]

intall = 2 (Integrate[i1, {x, 0, 1}, {y, -1, 0}] +
Integrate[i2, {x, 0, 1}, {y, 0, x}] +
Integrate[i3, {x, 0, 1}, {y, x, 1}]) // FullSimplify

(*   10/21 (7 (-1 + Sqrt[5]) \[Pi] - (
7 (-5 + Sqrt[5]) Sqrt[\[Pi]] Gamma[1/5])/(2^(1/5) Gamma[7/10]) -
70 Hypergeometric2F1[1/10, 7/10, 11/10, -1] +
10 Hypergeometric2F1[7/10, 7/10, 17/10, -1])   *)

intall // N

(*   35.2503   *)


it has to do with the order. If you integrate w.r.t. x first, then integrate the result w.r.t y then you get positive number.

If you integrate w.r.t. y first, then integrate the result w.r.t x then you get the imaginary number you showed.

ClearAll[x, y];
int = 1/(Abs[x]^(2/5) Abs[y]^(3/10) Abs[x - y]^(7/10));
N@Integrate[int, {y, -1, 1}, {x, -1, 1}]


N@Integrate[int, {x, -1, 1}, {y, -1, 1}]


Notice that Mathematica takes the integration variable in reverse order from the order they are listed.

The first one example above integrates over x first, then over y second, even though y appears first. And the second one integrates over y first, then over x second, even though x appears first. Which is what you had.

Now, why one gives positive real number and the second complex number? I do not know now. This needs more analysis. The Order of integration can make a difference of course (depending on the integrand).

But I do not know why one order gives complex value while the other order does not in this example.

• NIntegrate[int, {y, -1, 1}, {x, -1, 1}, Exclusions -> {x == y, x == 0, y == 0}, Method -> "AdaptiveMonteCarlo", AccuracyGoal -> 2, PrecisionGoal -> 2] produces 32.9916 and NIntegrate[int, {y, -1, 1}, {x, -1, 1}, Exclusions -> {x == y, x == 0, y == 0}, AccuracyGoal -> 3, PrecisionGoal -> 3] results in 35.2476 Oct 13 at 11:02
• Result of NIntegrate doesn't depend on the order of integration: NIntegrate[int, {y, -1, 1}, {x, -1, 1}, Exclusions -> {x == y, x == 0, y == 0}, AccuracyGoal -> 3, PrecisionGoal -> 3] and NIntegrate[int, {x, -1, 1}, {y, -1, 1}, Exclusions -> {x == y, x == 0, y == 0}, AccuracyGoal -> 3, PrecisionGoal -> 3] both give the same result 35.2476 Oct 13 at 11:18