A discrepancy between the analytical and numerical integration results

Question

Let us integrate the expression:

I*y*(1/(x - I*y)^(5/2) - 1/(x + I*y)^(5/2))

over a rectangle in the (x,y) plane, where x varies from -10 to 10, while y varies from 0 to 10. This integral can be solved analytically:

I*Integrate[
  y*(1/(x - I*y)^(5/2) - 1/(x + I*y)^(5/2)), {x, -10, 10}, {y, 0, 10}]
% // N

(* -(8/3) Sqrt[5] (2 Sqrt[2] + Sqrt[-1 + 5 Sqrt[2]] - \[Pi] - ArcCosh[3]) *)

(*  -2.31 *)

The last line here is the numeric value of the integral that I will use for comparisons below.

Now let us try to solve it numerically with several methods

 Quiet[Chop[
    I*NIntegrate[
      y *(1/(x - I y)^(5/2) - 1/(x + I y)^(5/2)), {x, -10, 10}, {y, 0,
        10}, Method -> #]]] & /@ {"LocalAdaptive", \
{"EvenOddSubdivision", Method -> "LocalAdaptive"}, 
  "AdaptiveMonteCarlo", "QuasiMonteCarlo", 
  "MonteCarlo", {"EvenOddSubdivision", 
   Method -> "AdaptiveMonteCarlo"}, {"EvenOddSubdivision", 
   Method -> "DuffyCoordinates"}}

(* {2.16, 2.16, 2.12, 2.26, 1.59, 2.33, 2.18} *)

Some of these methods do not, but some give warnings. I quieted them just to focus on the essential.

What strikes the eye here is that while the result of the exact solution is negative, the numerical result is positive.

Why? Any ideas?

**Edit: **

One can go to the cylindrical coordinates:

expr = Simplify[
  TrigToExp[
   I*y*(1/(x - I*y)^(5/2) - 1/(x + I*y)^(5/2)) /. {x -> 
      r*Cos[\[CurlyPhi]], y -> r*Sin[\[CurlyPhi]]}], {r > 0, 
   0 < \[CurlyPhi] < \[Pi]}]

(* (E^(-((7 I \[CurlyPhi])/
  2)) (-1 + E^(2 I \[CurlyPhi])) (-1 + E^(5 I \[CurlyPhi])))/(2 r^(
 3/2))  *)

and then integrate. In this case let us integrate it over the upper half-disk with the radius R=10:

Integrate[expr2*r, {r, 0, 10}, {\[CurlyPhi], 0, \[Pi]}] // N

(*   2.41   *)

and numerically

Quiet[NIntegrate[expr2*r, {r, 0, 10}, {\[CurlyPhi], 0, \[Pi]}, 
    Method -> #]] & /@ {"LocalAdaptive", {"EvenOddSubdivision", 
   Method -> "LocalAdaptive"}, "AdaptiveMonteCarlo", 
  "QuasiMonteCarlo", 
  "MonteCarlo", {"EvenOddSubdivision", 
   Method -> "AdaptiveMonteCarlo"}, {"EvenOddSubdivision", 
   Method -> "DuffyCoordinates"}}



 (*    {2.41, 2.41, 2.43, 2.42, 2.28, 2.38, 2.41}   *)

In this case they have the same sign.

Answer 1

Try to stay above y=0 by adding infinitesimally small positive value eps to lower integration limit for y and then take limit:

Limit[I*Integrate[y*(1/(x - I*y)^(5/2) - 1/(x + I*y)^(5/2)), {x, -10, 10}, {y, eps, 10}, Assumptions -> eps > 0], eps -> 0] // FullSimplify

(* (-4*Sqrt[10]*(-4 + Sqrt[-2 + 10*Sqrt[2]]))/3 *)

Answer 2

Integrate is wrong. Do ComplexExpand to get the right integral.

ceRe = FullSimplify[
   ComplexExpand[Re[I*y*(1/(x - I y)^(5/2) - 1/(x + I y)^(5/2))], 
   TargetFunctions -> {Re, Im}], y >= 0 && x \[Element] Reals]

NIntegrate[ceRe, {x, -10, 10}, {y, 0, 10}, WorkingPrecision -> 100, 
  MaxRecursion -> 50, AccuracyGoal -> 6, PrecisionGoal -> 6, 
  Exclusions -> {{0, 0}}]

(*   2.17329834703426526174107911927816608629584584526685923085449633834343\
       9718373530841484083310764572423   *)

cetr = TrigToExp[ceRe]

(*   (I y)/(2 ((x - I y)/Sqrt[x^2 + y^2])^(5/2) (x^2 + y^2)^(5/4)) - (
     I y ((x - I y)/Sqrt[x^2 + y^2])^(5/2))/(2 (x^2 + y^2)^(5/4)) - (
     I y)/(2 ((x + I y)/Sqrt[x^2 + y^2])^(5/2) (x^2 + y^2)^(5/4)) + (
     I y ((x + I y)/Sqrt[x^2 + y^2])^(5/2))/(2 (x^2 + y^2)^(5/4))   *)


Integrate[cetr, {x, -10, 10}, {y, 0, 10}]
   (N[#1, 20] &)[%]

(*   (11/606 + (3 I)/
 202) ((112 + 92 I) Sqrt[-10 + 10 I] + (22 - 18 I) Sqrt[-5 - 
 5 I] - (112 + 128 I) Sqrt[1 - 10 I] - (66 - 54 I) Sqrt[
 5 - 5 I] + (176 - 144 I) Sqrt[10] + (38 + 24 I) Sqrt[
 10 - 10 I] + (10 - 100 I) Sqrt[10 + 10 I] - (40 + 4 I) Sqrt[
 170 - 70 I] + 
 16 Sqrt[1105 + 262 I] - (44 - 36 I) Sqrt[
 5] \[Pi] - (40 + 4 I) Sqrt[10]
 ArcCoth[2 Sqrt[50/101 - (5 I)/101]] + (6 - 60 I) Sqrt[10]
 ArcCoth[Sqrt[1 + I]] - (88 - 72 I) Sqrt[5]
 ArcCoth[Sqrt[2]] - (6 - 60 I) Sqrt[10]
 ArcCoth[2/Sqrt[2 - I/5]] + (40 + 4 I) Sqrt[10]
 ArcCoth[(11/101 + (9 I)/101) Sqrt[10 - 100 I]] + (16 - 
  160 I) Sqrt[5] ArcTan[(1 + I)/Sqrt[2]] + (16 - 160 I) Sqrt[5]
 ArcTanh[(1 + I)/Sqrt[2]] + (2 - 20 I) Sqrt[10]
 Log[1 - 1/2 (-1 + I)^(3/2)] + (18 + 22 I) (-1)^(1/4) Sqrt[5]
 Log[1 + 1/2 (-1 + I)^(3/2)] + (2 - 20 I) Sqrt[10]
 Log[1 - 1/Sqrt[1 + I]] - (2 - 20 I) Sqrt[10]
 Log[1 + 1/Sqrt[1 + I]] - (3 - 30 I) Sqrt[10]
 Log[((1 + 2 I) + 2 Sqrt[-1 + I]) (3 - 2 Sqrt[2])] + (3 - 
  30 I) Sqrt[10] Log[2 - Sqrt[2]] - (3 - 30 I) Sqrt[10]
 Log[2 + Sqrt[2]] + (1 - 10 I) Sqrt[10]
 Log[-Sqrt[1 - I] + Sqrt[2]] - (1 - 10 I) Sqrt[10]
 Log[Sqrt[1 - I] + Sqrt[2]] - (10 + I) Sqrt[-5 - 5 I] Sqrt[1 - I]
 Log[-Sqrt[1 + I] + Sqrt[2]] + (10 + I) Sqrt[-5 - 5 I] Sqrt[1 - I]
 Log[Sqrt[1 + I] + Sqrt[2]] - (3 - 30 I) Sqrt[10]
 Log[(-(1/2) + I/2) Sqrt[1 + 10 I] + Sqrt[10]] + (3 - 30 I) Sqrt[
10] Log[(1/2 - I/2) Sqrt[1 + 10 I] + Sqrt[10]] - (1 - 10 I) Sqrt[
10] Log[2 Sqrt[5] - Sqrt[10 - I]] + (1 - I) Sqrt[100 - 495 I]
 Log[2 Sqrt[5] - Sqrt[10 - I]] + (1 - 10 I) Sqrt[10]
 Log[2 Sqrt[5] + Sqrt[10 - I]] - (1 - I) Sqrt[100 - 495 I]
 Log[2 Sqrt[5] + Sqrt[10 - I]])   *)

FullSimplify[%] gives

(*    8/3 (2 Sqrt[10] - Sqrt[-5 + 25 Sqrt[2]])    *)

(*  2.1732996388253877021 + 0.*10^-20 I  *)

Imaginary Part gives zero.

Answer 3

If you reverse the order of integration, you get the correct result:

I*Integrate[y*(1/(x - I*y)^(5/2) - 1/(x + I*y)^(5/2)), {y, 0, 10}, {x, -10, 10}]
N[%]
(*
8/3 Sqrt[5 (7 + 5 Sqrt[2] - 4 Sqrt[-2 + 10 Sqrt[2]])]
2.1733
*)

The difference probably has to do with the handling of the singularity at {0, 0}, such as described in What exactly does GenerateConditions do?

The difficulty is probably also connected with the branch cut for Power.