0
$\begingroup$

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
$\endgroup$

2 Answers 2

2
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ 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? $\endgroup$
    – yarchik
    Aug 17, 2021 at 12:31
  • $\begingroup$ Original post evaluates in Mathematica v12.2 0.250473 . Don't know which one is right. $\endgroup$ Aug 17, 2021 at 12:44
2
$\begingroup$

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}, 
 Method -> "LocalAdaptive"];

{f[0, 0, 0], f[1, 1, 0], f[1/2, 1/3, 1/9]}

(*   {0.340856, 0.0999031, 0.2748}   *)

pl1 = Plot[f[x, 0, 0], {x, -2, 2}, PlotPoints -> 15, 
  MaxRecursion -> 0]

Edit Now LessEqual in pw definition.

$\endgroup$
2
  • $\begingroup$ I tried to run your code in V 12.3, but I get the message "Unable to prove that integration limits {1,1/4 (-1-Sqrt[5]),...". Can you provide relevant assumptions? $\endgroup$ Aug 17, 2021 at 13:00
  • $\begingroup$ @AntonAntonov ,Does it work, if you Integrate ... over {t1, t2} ∈ region ? What do you mean with "relevant assumptions"? Maybe Element[{q1,q2,x,y,z} ,Reals]. In version 8.0 it works as shown. $\endgroup$
    – Akku14
    Aug 17, 2021 at 13:33

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.