# Nested NIntegrate over a region: Numerical Integration not converging

I need to computed two double nested numerical integrals, of which one is defined over a specific region (e.g. a pentagon). I've tried to use a single NIntegrate with four variables, but the result is highly oscillatory since it gives very different results for close evaluation points. How can I improve the calculation?

region = Polygon[Table[{Cos[2 Pi k/5], Sin[2 Pi k/5]}, {k, 5}]];

f[x_?NumericQ, y_?NumericQ, z_?NumericQ] := f[x, y, z] = (1/(2*Pi)^2)*NIntegrate[(1/(1 + q1^2 + q2^2))*Re[Exp[I*(z*Sqrt[1 + q1^2 + q2^2] + x*q1 + y*q2)]*Exp[I*(t1*q1 + t2*q2)]], {q1, -Infinity, +Infinity}, {q2, -Infinity, +Infinity}, {t1, t2} ∈ region];


Some of my (wrong) outputs are the following:

f[0, 0, 0]

-169.042 + 0. I

f[1, 1, 0]

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 4.34035 +0. I and 101.41637259371878 for the integral and error estimates.

0.109942 + 0. I


Just add Chop in your function:

f[x_?NumericQ, y_?NumericQ, z_?NumericQ] := (1/(2*Pi)^2)*
NIntegrate[(1/(1 + q1^2 + q2^2))*
Re[Exp[I*(z*Sqrt[1 + q1^2 + q2^2] + x*q1 + y*q2)]*
Exp[I*(t1*q1 +
t2*q2)]], {q1, -Infinity, +Infinity}, {q2, -Infinity, \
+Infinity}, {t1, t2} \[Element] region ] // Chop;

f[0,0,0] //AbsoluteTiming
(*{20.6637, 0.250473}*)


But evaluation takes some time ~20s.

• With your code I am getting f[0,0,0]==0.208878, in the original post f[0, 0, 0]==-169.042 + 0. I, @Akku14 reports f[0,0,0]==0.340856. Which one is right? Aug 17, 2021 at 12:31
• Original post evaluates in Mathematica v12.2 0.250473 . Don't know which one is right. Aug 17, 2021 at 12:44

You can do the t1,t2 integration anylytically, which make things much faster and stable.

region = Polygon[Table[{Cos[2 Pi k/5], Sin[2 Pi k/5]}, {k, 5}]]


Since i am working with MMA version 8.0, i can not use {t1, t2} ∈ region, but had to construct a piecewise function pw.

{{x1, y1}, {x2, y2}, {x3, y3}} = {{1/4 (-1 - Sqrt[5]), Sqrt[
5/8 - Sqrt[5]/8]}, {1/4 (-1 + Sqrt[5]), Sqrt[
5/8 + Sqrt[5]/8]}, {1, 0}};

pw[t1_] =
Piecewise[{{y1 + (t1 - x1) *(y2 - y1)/(x2 - x1),
x1 <= t1 < x2}, {y2 + (t1 - x2)*(y3 - y2)/(x3 - x2),
x2 <= t1 <= x3}}] // Simplify

Plot[-pw[t1], {t1, -1, 1}]

Graphics[region, Frame -> True]

integrand = (1/(2*Pi)^2)*(1/(1 + q1^2 + q2^2))*
Re[Exp[I*(z*Sqrt[1 + q1^2 + q2^2] + x*q1 + y*q2)]*
Exp[I*(t1*q1 + t2*q2)]] //
ComplexExpand[#, TargetFunctions -> {Re, Im}] & // Simplify

(*   Cos[q1 (t1 + x) + q2 (t2 + y) +
Sqrt[1 + q1^2 + q2^2] z]/(4 \[Pi]^2 (1 + q1^2 + q2^2))   *)

int[q1_, q2_, x_, y_, z_] =
Integrate[Cos[q1 (t1 + x) + q2 (t2 + y) + Sqrt[1 + q1^2 + q2^2] z]/(
4 \[Pi]^2 (1 + q1^2 + q2^2)), {t1, 1/4 (-1 - Sqrt[5]),
1}, {t2, -pw[t1], pw[t1]}] // Simplify
(*   Takes 2 minutes   *)

f[x_?NumericQ, y_?NumericQ, z_?NumericQ] :=
NIntegrate[
int[q1, q2, x, y,
z], {q1, -Infinity, +Infinity}, {q2, -Infinity, +Infinity},
`