1

I am using NIntegrate to evaluate the integral of the real part of some function.

nbs1 = {{0, 1}, {-1, 1}, {0, 0}};
A1 = 0.6049; l = 3 A1 + 10^-3;
gap[k_] := A1*Sum[Exp[-2 π I k . nb], {nb, nbs1}];
en[k_] := Sqrt[l^2 - Abs[gap[k]]^2]; 
fg[k_] := gap[k]/(l + en[k]);

NIntegrate[Re[fg[{kx, ky}] Exp[2 π I (kx + ky)]],
 {kx, -0.5, 0.5}, {ky, -0.5, 0.5}]

However, the result is 0.014725 + 0. I, still carrying a complex part although it is identically zero, and I also get the message "Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small".

But the integrand is explicitly real (by the function Re; this contrasts some old posts like this and this one), and I am using NIntegrate instead of Integrate. So I wonder what produces this + 0. I (how Mathematica calculate the integral under the hood), and how I can get rid of it during the integration and possibly improve the evaluation speed.

Zhengyuan Yue
  • 179
  • 5
  • You can just Chop the integral result to get rid of the + 0. I. You can also ComplexExpand the integrand, although this does make the integration take longer:

    int = Re[fg[{kx, ky}] Exp[2 \[Pi] I (kx + ky)]] // ComplexExpand; NIntegrate[int, {kx, -0.5, 0.5}, {ky, -0.5, 0.5}] (*0.014725*)

    The NIntegrate warnings still persist however

     – ydd Sep 03 '23 at 14:43
  • 1
    I am aware of these methods; however, I want to know more on how Mathematica do my integral under the hood. – Zhengyuan Yue Sep 03 '23 at 14:46
  • Ah ok. For what it's worth, explicitly stating the Method as Method -> "LocalAdaptive" seems to remove this problem, while Method -> "GlobalAdaptive" has the +0. I part and the warning messages. – ydd Sep 03 '23 at 14:52
  • Re[NIntegrate[ Re[fg[{kx, ky}] Exp[2 \[Pi] I (kx + ky)]], {kx, -0.5, 0.5}, {ky, -0.5, 0.5}]] ? – Mariusz Iwaniuk Sep 03 '23 at 15:28
  • Use A1 = 0.6049 // Rationalize; and NIntegrate[Re[fg[{kx, ky}] Exp[2 \[Pi] I (kx + ky)]], {kx, -1/2, 1/2}, {ky, -1/2, 1/2}, WorkingPrecision -> 15, MinRecursion -> 10, MaxRecursion -> 20] – Bob Hanlon Sep 03 '23 at 15:57

1 Answers1

3

Overview

  1. The problem arises in the compiled function created for some sort of singularity processing that NIntegrate[] decided was necessary.

  2. Is it a bug? Ask WRI. Here are two things to consider. (a) The singularity processing is not necessary. (b) The compiled function created has a return type Complex even though the function is Real. The function being compiled is derived from the integrand and is slightly more complicated (a linear combination of Re[] expressions). Compile[] seems to gives up on the analysis of this function and defaults to the safer type Complex for the result.

  3. The comments show some ways to avoid a complex output. A direct fix is to turn off the singularity handler:

NIntegrate[
 Re[fg[{kx, ky}] Exp[2 \[Pi] I (kx + ky)]],
  {kx, -0.5, 0.5}, {ky, -0.5, 0.5},
 MinRecursion -> 1, (* stops the GlobalAdaptive slow-con warning *)
 Method -> {"GlobalAdaptive", "SingularityHandler" -> None}]
(*  0.014725  *)

Analysis

Here's how to analyze the integration. IntegrationMonitor returns information about the integral over each subregion.

NIntegrate[
 Re[fg[{kx, ky}] Exp[2 \[Pi] I (kx + ky)]],
 {kx, -0.5, 0.5}, {ky, -0.5, 0.5},
 IntegrationMonitor :> ((foo = #) &)]

Pick out the subregions over which the integral is complex:

Position[Through[foo["Integral"]], _Complex]
cplx = Extract[foo, %];
(*
{{137}, {138}, {139}, {140}, {141}, {142}, {146}, {147}, {148}, \
{149}, {150}, {155}, {156}, {327}, {328}, {333}, {336}, {337}, {338}, \
{339}, {342}, {344}, {346}, {347}, {348}, {349}}
*)

Each subregion has an Experimental`NumericalFunction[] used to evaluate the integral, and the numerical function has a CompiledFunction (in the OP's case). Unsurprisingly, the ones related to Real integral values have Real integrand values, and those with Complex integral values have Complex integrand values.

nfR = foo[[1]]["NumericalFunction"];
cfR = nfR@"CompiledFunction";

nfC = cplx[[1]]["NumericalFunction"];
cfC = nfC@"CompiledFunction";


nfR[-0.375, -0.375]
cfR[-0.375, -0.375]
(*
-0.125894
-0.125894
*)


nfC[0.96875`, 1.]
cfC[0.96875`, 1.]
(*
0.00102451 + 0. I
0.00102451 + 0. I
*)

We can check that the return type agrees with the results (a redundant step).

Needs@"CompiledFunctionTools`"

Cases[ToCompiledProcedure[cfR], _CompiledResult, Infinity]
(*  {CompiledResult[Register[Real, 2]]}  *)


Cases[ToCompiledProcedure[cfC], _CompiledResult, Infinity]
(*  {CompiledResult[Register[Complex, 17]]}  *)

One can examine the end of the output of CompilePrint@cfR and CompilePrint@cfC to see that both apply Re[] etc. The output is rather long, and I will omit it. One can also see whether the result is a real or complex register.

We can test the hypothesis that it is Compile that fails to parse the Real value of the function:

testC = Compile[{{kx, _Complex}, {ky, _Complex}},
  Evaluate@nfC@"FunctionExpression"]
testC[0.96875`, 1.]
(*  0.00102451 + 0. I  *)

Here's the form of the function in nfC[]:

Shallow[nfC@"FunctionExpression", 5]
(*  kx (0.00236289 Re[<<1>>] + 0.00236289 Re[<<1>>])  *)

Clearly, one wants the time Compile[] spends analyzing code to be very short, especially in auto-compiled code. It's hard to argue strongly that it is a bug. Normally, if the user is calling Compile[] there are workarounds. I don't know of any hooks in NIntegrate[] to control Compile[].

Michael E2
  • 235,386
  • 17
  • 334
  • 747