4
$\begingroup$

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.

$\endgroup$
12
  • 1
    $\begingroup$ Maybe a different branch on the Sqrt (i.e., the ^5/2)? $\endgroup$
    – bill s
    Commented Feb 27, 2018 at 15:25
  • $\begingroup$ @ bill s Have you an idea of how to fix the branch? $\endgroup$ Commented Feb 27, 2018 at 15:30
  • $\begingroup$ Maple is giving me 2.173299654 $\endgroup$
    – zhk
    Commented Feb 27, 2018 at 15:30
  • $\begingroup$ Even stranger: you can check that the integrand is equal to $f(x,y)$ with f[x_, y_] = y/(x^2 + y^2)^(5/2)*((x + I*y)^(5/2) - (x - I*y)^(5/2)). However, I*Integrate[f[x, y], {x, -10, 10}, {y, 0, 10}] returns -(8/3) Sqrt[-5 + 25 Sqrt[2]] that is approx. -14.69. Since there is a singularity in $(0,0)$ I guess it might be related to the chosen branches. (That's probably not a problem here but note that the integrals to note commute because of the singularity). $\endgroup$
    – anderstood
    Commented Feb 27, 2018 at 15:31
  • $\begingroup$ @bills Are there branches even for surface integral, btw? Does Mathematica always uses Stokes theorem to solve a contour integral? $\endgroup$
    – anderstood
    Commented Feb 27, 2018 at 15:32

3 Answers 3

6
$\begingroup$

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 *)

$\endgroup$
4
$\begingroup$

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.

$\endgroup$
5
  • 1
    $\begingroup$ Trying to verify your code the last Integrate evaluates to -14.6922???(MMA 11.0.1 Windows64) $\endgroup$ Commented Feb 27, 2018 at 17:00
  • $\begingroup$ To @Ulrich Neumann. I am working with MMA Version 8.0. Look at the full output of Integrate I added. $\endgroup$
    – Akku14
    Commented Feb 27, 2018 at 17:05
  • 1
    $\begingroup$ The result of Integrate (MMA 11.0.1) is significant shorter(and obviosly wrong): (* -(8/3) Sqrt[-5 + 25 Sqrt[2]]*) $\endgroup$ Commented Feb 27, 2018 at 17:11
  • $\begingroup$ FullSimplified it is also much shorter 8/3 (2 Sqrt[10] - Sqrt[-5 + 25 Sqrt[2]]) . $\endgroup$
    – Akku14
    Commented Feb 27, 2018 at 17:18
  • $\begingroup$ @UlrichNeumann Same value as in my comment above. $\endgroup$
    – anderstood
    Commented Feb 27, 2018 at 18:21
3
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ What allows you to reverse the order of integration? See for instance with Integrate[(x^2 - y^2)/(x^2 + y^2)^2, {x, 1, Infinity}, {y, 1, Infinity}] (taken from wikipedia). $\endgroup$
    – anderstood
    Commented Feb 28, 2018 at 16:30
  • 2
    $\begingroup$ @anderstood The region is compact and outside a (half) disk of radius $\epsilon$ at the origin the function is bounded; and the integral over the half disk is $16 \sqrt{\epsilon}\,/\,21$, and so vanishes as $\epsilon \rightarrow 0$. Hence, with the aid of Fubini's theorem, one can show the order does not matter. (In polar coordinates, the integrand is rather easy to understand.) $\endgroup$
    – Michael E2
    Commented Feb 28, 2018 at 17:23
  • $\begingroup$ Thank you for the explanation. $\endgroup$
    – anderstood
    Commented Mar 1, 2018 at 17:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.