# Error with NIntegrate of a two variables integration

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.

• "I include some poles in the integration range" - what do you mean by that? Can you include the code you used? Mar 18, 2021 at 14:36
• I don't see any singularities on the domain of the integration since the denominator equals zero if y == (-2 - I) - x || y == -I - x || y == I - x || y == (2 + I) - x || y == (-2 - I) + x || y == -I + x || y == I + x || y == (2 + I) + x. I don't know good numerical methods for multple improper integras. Try to cut the tails, replacing Infinity by 100. Mar 18, 2021 at 16:40
• Thanks for your comments. "I include some poles in the integration range" just mean I change the integration range as for example, {y, 0,2+x,-Infinity}, here I ignore the imaginary part. And this is an integration along real axis, i don't know why the poles with imaginary parts maybe hurt the integration. I trust the integration will be converged due to the power counting of x, as @user64494 said. Mar 19, 2021 at 1:45

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.

• If I correctly count, the total degree of the numerator equals 4 and the total degree of the denominator equals 12. In view of it I believe that the double integral under consideration converges. Mar 18, 2021 at 17:37
• Sounds so for sure. Would be nice to understand why my code is not exhibiting that. Mar 18, 2021 at 17:52
• Look at Plot3D[{Im[((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, 100}, {x, -100, 0}, PlotRange -> All, PlotPoints -> 100]. Maybe, the ravine along y==-x is badly handled. Mar 18, 2021 at 18:10
• +1. No comments. Mar 19, 2021 at 6:23
• I get an enormous speed up of your code when: removing everywhere MaxRecursion -> 25, WorkingPrecision -> 30; changing the integration method to Method -> "LocalAdaptive". Then the result I get are (*-2.6879624115154894*) and (*-3.134223337637785*) Mar 19, 2021 at 9:20

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!

• I try to use your methods to investigate the new function myf(see update), and the real part times r will also tend to a constant as r increases to infinity. But this new function will be integrated as 1/pi. Mar 19, 2021 at 12:58
• @Jiang With myf[] there is only a real part of r ReIm[i\[CurlyPhi][r] ] , integral doesn't converge ! By the way this case can be solved by Mathematica and gives a message ...does not converge on {0,\[Infinity]} Mar 19, 2021 at 13:06
• I get your point, actually the integration of myf is dependent of order of integral, the result I post is just integrating over x and then y, but if we change the order, the integration will not be converged. I think the first integration maybe have the same question so the MMA can not give the numerical results. Thanks for your answer! Mar 19, 2021 at 14:31
• @Jiang That's the point! You're welcome. Mar 19, 2021 at 14:35
• +1. No comments. Mar 19, 2021 at 15:58

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

• Crash occurs after about 22 minutes on 12.1 (Linux 64-bit) after using about 25GB of RAM. Mar 20, 2021 at 14:43
• I had to leave my computer (Mac) for several minutes, and when I came back there was the running-out-of-memory dialog and several (but not all) apps had crashed, including Mma (kernel & FE). MemoryInUse[] had gotten up to 12GB when I left. Mar 20, 2021 at 16:00
• @eyorble and MichaelE2 -- Thank you! I will include your reports in my bug report to WRI. Mar 20, 2021 at 17:00