Error with NIntegrate of a two variables integration

Question

This is mathematica code for my integration,

NIntegrate[((8 - 8 I) x y (1 + x^2 + y^2) ((3 + 4 I) + (2 + 4 I) x^2 -
     x^4 - 2 ((1 + 2 I) + x^2) y^2 + 
    3 y^4))/(\[Pi] (-((-2 - I) + x)^2 + y^2) (-((2 + I) + x)^2 + 
    y^2) ((1 + x^2)^2 - 2 (-1 + x^2) y^2 + y^4)^2), {y, 0, 
  Infinity}, {x, -Infinity, 0}]

where the denominator has some poles which have imaginary part. The code gives an error result,

-94.2449 - 98.463 I

and the error messages are

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 18 recursive bisections in y near {y,x} = {3008.02,-3007.26}. NIntegrate obtained -94.2449-98.463 I and 26.019110660515846` for the integral and error estimates.

Then I try increase MaxRecursion

NIntegrate[((8 - 8 I) x y (1 + x^2 + y^2) ((3 + 4 I) + (2 + 4 I) x^2 -
     x^4 - 2 ((1 + 2 I) + x^2) y^2 + 
    3 y^4))/(\[Pi] (-((-2 - I) + x)^2 + y^2) (-((2 + I) + x)^2 + 
    y^2) ((1 + x^2)^2 - 2 (-1 + x^2) y^2 + y^4)^2), {y, 0, 
  Infinity}, {x, -Infinity, 0}, MaxRecursion -> 40]

But it doesn't improve.

66.4382 + 68.1894 I

and the messages are

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 66.4382 +68.1894 I and 204.94038847839664` for the integral and error estimates.

Finally, I include some poles in the integration range, but it does not help. So my questions are

1. how to get the correct answer of this integration?
2. Are there any general methods to solve Multiple integrals which have poles with Mathematica?

Update

I think the integration will be inefficient near the poles. Here I try another example,

myf[x_, y_] := (
 8 x y (1 + x^2 + y^2) (-(1 + x^2)^2 - 2 (-1 + x^2) y^2 + 
    3 y^4))/(\[Pi] ((1 + x^2)^2 - 2 (-1 + x^2) y^2 + y^4)^3)
NIntegrate[myf[x, y], {x, -Infinity, 0}, {y, 0, Infinity}]

It also gives the bad result

out:-105.399

But this function can be integrated analytically

Integrate[myf[x, y], {x, -Infinity, 0}]

out:ConditionalExpression[(
 2 y)/(\[Pi] (1 + y^2)^2), -1 < Im[y] < 1 || Im[y] > 1 || Im[y] < -1]

and then

Integrate[(2 y)/(\[Pi] (1 + y^2)^2), {y, 0, Infinity}]

out:1/\[Pi]

From the denominator the poles can be solved as

{{x -> -I - y}, {x -> -I - y}, {x -> -I - y}, {x -> I - y}, {x -> 
   I - y}, {x -> 
   I - y}, {x -> -I + y}, {x -> -I + y}, {x -> -I + y}, {x -> 
   I + y}, {x -> I + y}, {x -> I + y}}

I guess we can get some ideas from firstly find a numerical solution of this example.

Update

The integration of myf will be dependent of order of integral. This maybe the reasons why MMA can't give the correct answer.

Answer 1

See update below as per user64494

Also, see update to integration method suggested by Daniele below

Do you know beforehand if it converges? I can create an integral function innerI(x) of the inner integral from say y from 0 to 2000 as a function of x, then integrate that integral from say x from 0 to -400 and the absolute value doesn't appear to be converging. You may wish to further study my code for larger ranges and perhaps the trend may settle down:

myf[x_, y_] := ((8 - 8 I) x y (1 + x^2 + 
       y^2) ((3 + 4 I) + (2 + 4 I) x^2 - x^4 - 
       2 ((1 + 2 I) + x^2) y^2 + 3 y^4))/(\[Pi] (-((-2 - I) + x)^2 + 
       y^2) (-((2 + I) + x)^2 + 
       y^2) ((1 + x^2)^2 - 2 (-1 + x^2) y^2 + y^4)^2);
(*
 create integral function of inner integral for y from 0 to 2000
*)
innerI[theX_?NumericQ] := 
  NIntegrate[myf[theX, y], {y, 0, 2000}, MaxRecursion -> 20];
(*
  generate table of integral of innerI[x] for x from 0 to -400
*)
theTable = Table[
   {xMin, Abs@NIntegrate[innerI[x], {x, xMin, 0}]},
   {xMin, 0, -400, -10}];
ListPlot[theTable, Joined -> True]

Update:

Based on user64494's analysis of the integral converging, I'd like to try again. This time, I integrate over either sides of the ravines as I think an argument could be made the function decays rapidly away from them: First define function:

myf[x_, y_] := ((8 - 8 I) x y (1 + x^2 + 
       y^2) ((3 + 4 I) + (2 + 4 I) x^2 - x^4 - 
       2 ((1 + 2 I) + x^2) y^2 + 3 y^4))/(\[Pi] (-((-2 - I) + x)^2 + 
       y^2) (-((2 + I) + x)^2 + 
       y^2) ((1 + x^2)^2 - 2 (-1 + x^2) y^2 + y^4)^2);

(*
  set integration -delta+ravine to ravine+50 as an approximation
*)
delta = 50;
(*
 first integrate -50<x<0, 0<y50
*)
myf1[currentY_?NumericQ] := 
 NIntegrate[myf[x, currentY], {x, -50, 0}, MaxRecursion -> 25, 
  WorkingPrecision -> 30]
int1 = NIntegrate[myf1[y], {y, 0, 50}, MaxRecursion -> 25, 
   WorkingPrecision -> 30];
(*
 next, set up inner integral to integrate along the primary trough \
y=-x from -(y+50) to -(y-50) 
*)
myf2[currentY_?NumericQ] := 
  NIntegrate[
   myf[x, currentY], {x, -(currentY + delta), -(currentY - delta)}, 
   MaxRecursion -> 25, WorkingPrecision -> 30];
(*
 generate some points y=50 to 500
*)
table2 = Table[
   {maxY, 
    NIntegrate[myf2[y], {y, 50, maxY}, MaxRecursion -> 25, 
     WorkingPrecision -> 30]},
   {maxY, 50, 500, 50}];

(*
 separate real and imag parts
*)
reTrend = {#[[1]], Re[int1 + #[[2]]]} & /@ table2;
imTrend = {#[[1]], Im[int1 + #[[2]]]} & /@ table2;
(*
 get final results at y=500 as approximate value of integral
*)
Last@reTrend[[All, 2]]
Last@imTrend[[All, 2]]
(*
 plot trend
*)
ListPlot[{reTrend, imTrend}, PlotStyle -> {Red, Blue}, Joined -> True]

-2.68796241386416885973535164513

-3.13422333638947755352380311449

Basically, I'm suggesting the value of the integral is near $-2.6879-3.134i$ and would be curious what other interested persons obtain.

Answer 2

If the integral converges it is allowed to change the coordinates {x,y}->{r,\[CurlyPhi]} to polar and to consider the eqivalent integral Integrate[r "integrand",{r,0,Infinity},{\[CurlyPhi],Pi/2,Pi}]!

integrand=  ((8 - 8 I) x y (1 + x^2 + y^2) ((3 + 4 I) + (2 + 4 I) x^2 -x^4 - 2 ((1 + 2 I) + x^2) y^2 +3 y^4))/(\[Pi] (-((-2 - I) + x)^2 + y^2) (-((2 + I) + x)^2 + y^2) ((1 + x^2)^2 - 2 (-1 + x^2) y^2 + y^4)^2) 
/. {x ->r Cos[\[CurlyPhi]], y -> r Sin[\[CurlyPhi]]} // FullSimplify

(*((16 - 16 I) r^2 (1 + r^2) ((6 + 8 I) + r^4 -4 r^2 ((-1 - 2 I) + r^2) Cos[2 \[CurlyPhi]] +r^4 Cos[4 \[CurlyPhi]]) Sin[2 \[CurlyPhi]])
/(\[Pi] ((-14 + 48 I) - (12 + 16 I) r^2 + r^4 + r^4 Cos[4 \[CurlyPhi]]) (2 + 4 r^2 + r^4 +r^4 Cos[4 \[CurlyPhi]])^2)*)

Now we split the integration to Integrate[r Integrate["integrand",{\[CurlyPhi],Pi/2,Pi}],{r,0,Infinity}] and look at the inner integral (over \[CurlyPhi])

i\[CurlyPhi][r_?NumericQ] :=NIntegrate[integrand, {\[CurlyPhi], Pi/2, Pi} , Method -> "LocalAdaptive"  ]

This result times r forms the integrand for the pending integration (over r) and looks like (real and imaginary part)

Plot[r ReIm[i\[CurlyPhi][r] ] // Evaluate, {r, 0, 100},PlotRange -> All, PlotStyle -> {Blue, Red}, AxesLabel -> {"r", "re[r],im[r]"}]  

Both curves tend to a constant value for r->Infinity , that means the integration Integrate[ u ReIm[i\[CurlyPhi][u] ],{u,0,r}] is increasing more and more !

That's why the integral OP asked for doesn't converge, I think!

Answer 3

( Not an answer, a crash proclaim and verification request. )

The code below crashes the Mathematica kernel on my computer, Version 12.2.0 for Mac OS X x86 (64-bit) (December 12, 2020).

Please, run it on your machines/versions and proclaim the results. E.g. on Windows 10, Mathematica Version 12.1, etc. (It takes 20-35 min to get the crash.)

AbsoluteTiming[
 NIntegrate[((8 - 8 I) x y (1 + x^2 + 
       y^2) ((3 + 4 I) + (2 + 4 I) x^2 - x^4 - 
       2 ((1 + 2 I) + x^2) y^2 + 3 y^4))/(\[Pi] (-((-2 - I) + x)^2 + 
       y^2) (-((2 + I) + x)^2 + 
       y^2) ((1 + x^2)^2 - 2 (-1 + x^2) y^2 + y^4)^2), {y, 0, 
   Infinity}, {x, -Infinity, 0}, MaxRecursion -> 100, 
  Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 100000}]
 ]

(*NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.*)