6
$\begingroup$

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.
$\endgroup$
1
  • 2
    $\begingroup$ The same issue in 12.3.1 on Windows 10. RealAbs instead of Abs does not help. $\endgroup$
    – user64494
    Oct 13 at 9:35
7
$\begingroup$

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   *)
$\endgroup$
6
$\begingroup$

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

Mathematica graphics

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

Mathematica graphics

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.

$\endgroup$
3
  • $\begingroup$ 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 $\endgroup$
    – user64494
    Oct 13 at 11:02
  • 3
    $\begingroup$ 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 $\endgroup$ Oct 13 at 11:18
  • 1
    $\begingroup$ @Nasser, please see my answer. $\endgroup$
    – Akku14
    Oct 13 at 14:38

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.