1
$\begingroup$

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.

Improve this question
$\endgroup$
5
  • $\begingroup$ 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 $\endgroup$ Commented Sep 3, 2023 at 14:43
  • 1
    $\begingroup$ I am aware of these methods; however, I want to know more on how Mathematica do my integral under the hood. $\endgroup$ Commented Sep 3, 2023 at 14:46
  • $\begingroup$ 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. $\endgroup$ Commented Sep 3, 2023 at 14:52
  • $\begingroup$ Re[NIntegrate[ Re[fg[{kx, ky}] Exp[2 \[Pi] I (kx + ky)]], {kx, -0.5, 0.5}, {ky, -0.5, 0.5}]] ? $\endgroup$ Commented Sep 3, 2023 at 15:28
  • $\begingroup$ 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] $\endgroup$ Commented Sep 3, 2023 at 15:57

1 Answer 1

Reset to default
3
$\begingroup$

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[].

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.