Fn[x_, y_, z_, r_, 
   br_] = -1536 Im[(
    6 Cos[z] + 
     2 Cos[x] (3 + (-7 + 3 Cos[y]) Cos[z]))/(478 + 96 I br + 
      32 br^2 + 96 r - 64 I br r - 32 r^2 + 11 Cos[2 x] - 
      264 Cos[y] - 48 I br Cos[y] - 48 r Cos[y] + 11 Cos[2 y] - 
      336 Cos[z] + 144 Cos[y] Cos[z] + 
      12 Cos[x] (-22 - 4 I br - 4 r + 9 Cos[y] + 12 Cos[z]))^2];  

I have the above function and I would like to perform numerical integration using NIntegrate with respect to (x,y,z) and then r with different upper limits scaled as R as follows

sg = ParallelTable[{R, 
      Fn[x, y, z, r, 
       0.01], {x, -\[Pi], \[Pi]}, {y, -\[Pi], \[Pi]}, {z, -\[Pi], \
\[Pi]}, {r, -5, R}, 
      Method -> {Automatic, "SymbolicProcessing" -> 0}]}, {R, -1, 1.5,
      0.05}] // AbsoluteTiming;

when I plot the results this is what I get

ListLinePlot[sg[[2]][[All, {2, 1}]], 
 PlotStyle -> {{Red, Thickness[0.01]}}, Frame -> True, Axes -> False, 
 FrameLabel -> {"F(R)", "R"}, 
 LabelStyle -> {FontFamily -> "Latin Modern Roman", Black, 
   FontSize -> 16}, PlotRange -> {Full, Full}, ImageSize -> 300, 
 AspectRatio -> 1]   

enter image description here

while the expected form of the curve must be like this (updated)
enter image description here
So, how can I improve the numerical integration to obtain the desired results?

  • $\begingroup$ 1. Your reference data has no labels on the x-axis. 2. The integrand function has singularities. $\endgroup$
    – yarchik
    Jul 26, 2020 at 13:52
  • $\begingroup$ @yarchik then what? I mean if it has singularities can't I use NIntegrate ? $\endgroup$ Jul 26, 2020 at 13:56
  • 1
    $\begingroup$ If an integrand has singularities then the result is generally undefined. There are many cases when a well-defined value can be assigned to the integral. In 1D one would compute a Cauchy principal value. For multidimensional integrals this option is not operative and one needs to more closely analyze the singularities. You can specify them with Exclusions options (reference.wolfram.com/language/tutorial/…) $\endgroup$
    – yarchik
    Jul 26, 2020 at 14:06
  • 1
    $\begingroup$ The integrand does have singularities, but I don't think you ever encounter them in the given integration region. The denominator has both a real and imaginary component, and both Reduce and FindInstance tell me that there is not point in the integration region where both components vanish simultaneously, so we should be fine. $\endgroup$
    – Hausdorff
    Jul 28, 2020 at 12:41
  • 1
    $\begingroup$ I tried some integration strategies for your previous question, and I got results similar to what you are expecting using the AdaptiveMonteCarlo method. I needed a huge number of points, though, and the accuracy was still pretty bad. The LocalAdaptive strategy also gave me results, but they disagreed with the Monte Carlo ones, and increasing MinRecursion something like 19 did not change anything. For performance I you also could compile the function, that way you can get a speedup of another order of magnitude or two. $\endgroup$
    – Hausdorff
    Jul 28, 2020 at 12:44

2 Answers 2


First of all, there do not seem to be any singularities in the integration region:

$Assumptions = And[-Pi<=x<=Pi,Pi<=y<=Pi,Pi<=z<=Pi,-5<=r<=2];

den = 478 + (96*I)*br + 32*br^2 + 96*r - (64*I)*br*r - 32*r^2 + 11*Cos[2*x] - 
      264*Cos[y] - (48*I)*br*Cos[y] - 48*r*Cos[y] + 11*Cos[2*y] - 336*Cos[z] + 
      144*Cos[y]*Cos[z] + 12*Cos[x]*(-22 - (4*I)*br - 4*r + 9*Cos[y] + 12*Cos[z]) /. {br->1/100};

denReIm = Through[{Re,Im}[den]] // FullSimplify;

Reduce[denReIm=={0,0} && $Assumptions, {x,y,z,r}]
FindInstance[denReIm=={0,0} && $Assumptions, {x,y,z,r}]

By compiling Fn, you can further improve the evaluation time by a factor of about 100

FnCompiled = Compile[{{x,_Real},{y,_Real},{z,_Real},{r,_Real},{br,_Real}},Evaluate@Fn[x,y,z,r,br]]

Fn@@RandomReal[{-Pi,Pi},{5}] // RepeatedTiming
FnCompiled@@RandomReal[{-Pi,Pi},{5}] // RepeatedTiming
{0.0000176, -0.00121113}
{2.428*10^-6, -0.0177212}

The integration seems to converge to something sensible using the AdaptiveMonteCarlo method:

FvMC[v_,points_] := NIntegrate[FnCompiled[x,y,z,r,1/100],{r,-5,v},{x,-Pi,Pi},{y,-Pi,Pi},{z,-Pi,Pi},Method->{"AdaptiveMonteCarlo",MaxPoints->points}];

data7 = ParallelTable[{v,FvMC[v,10^7]}, {v, -1, 1.5, 5/36}];
data8 = ParallelTable[{v,FvMC[v,10^8]}, {v, -1, 1.5, 5/34}];

    MapAt[Around[#,14]&, data7, {All,2}],
    MapAt[Around[#,4.2]&, data8, {All,2}]
},IntervalMarkers->"Bands",PlotLegends -> {"10^7 Points","10^8 Points"}]

enter image description here

The error bands are based on the error estimates Mathematica gave me in the NIntegrate::maxp errors. For both runs, the average time per point appears to be around $1.4\times10^{-5}$ seconds, which is quite a bit slower than what it should be, but I suppose that could be attributed to overhead.


With Akku14's suggestions in the comments, i.e splitting up the integration region into subintervals, as well as using the symmetries $(x\leftrightarrow -x)$, $(y\leftrightarrow -y)$, $(z\leftrightarrow -z)$ of the integrand, we can further improve the result. The integration method LocalAdaptive now also seems to be able to give results in a reasonable time frame, so I'll also include it below.

We have to modify the integration functions to be



In the LocalAdaptive case it is important to set MinRecursion to something larger than 3. While NIntegrate does not throw an error for lower values, both the result and runtime make a jump, and tend to remain stable from there on out (in the following I am always using MinRecursion->15). My guess would be that there are features in the integrand that are too small for NIntegrate to notice for small values, but I am not really sure what is happening.

Performing the integrations in steps of 0.05 using the LocalAdaptive method, as well as AdaptiveMonteCarlo with both 10^7 and 10^8 points per interval, I get the following result:

enter image description here

The Monte Carlo results appears to converge to that of the LocalAdaptive strategy. The difference increases for larger values of v, but that is to be expected, as we are adding all the errors of the previous steps.

The plot is normalized to the point {-1,0}, since AdaptiveMonteCarlo and LocalAdpative appear to disagree on the value of the integral from -5 or -1,

NIntegrate::maxp: The integral failed to converge after 100000100 integrand evaluations. NIntegrate obtained -5.3294 and 0.31752764339093775` for the integral and error estimates.

My guess would be that this is caused by features that are overlooked by NIntegrate, since the integrand is rather 'spikey'. Perhaps splitting up this integration region into smaller subintervals may help here as well.

Also note that even with the performance improvements the LocalAdaptive strategy is still very time consuming. The following graph shows the required time per interval of the results above: strong text

  • $\begingroup$ this looks great! Note that r runs from -5 (not -2) to v $\endgroup$ Jul 28, 2020 at 14:23
  • $\begingroup$ Sorry, I just fixed it $\endgroup$
    – Hausdorff
    Jul 28, 2020 at 16:44
  • 2
    $\begingroup$ Dear @Hausdorff, first, you can reduce time for integration by a factor of 8, if you regard Fn[x, y, z, r] == Fn[-x, y, z, r] == Fn[x, -y, z, r] == Fn[x, y, -z, r] == Fn[-x, -y, -z, r] . Second, do integration from r==-5 to r==-1 only once, store it and integrate then from -1 to -1+5/36, then from -1+5/36 to -1+2*5/36 and so on. Saves a lot of time. $\endgroup$
    – Akku14
    Jul 29, 2020 at 3:21
  • $\begingroup$ Thank you, you are completely right. I'll try to work both of your points in later today $\endgroup$
    – Hausdorff
    Jul 29, 2020 at 10:11

Based on the detailed answer by @Hausdorff and the hints from @Akku14, I found that it is more efficient to transform the space integral into ParallelSum and then perform NIntegrate at each point with respect to r. Now the whole process takes only 15 min for all steps and with acceptable accuracy as shown in the attach Fig. Note that the max. value is about 60 (more close to the real one shown in the question) in the Fig. below while it was almost doubled, about 120 in the last result of @Hausdorff.

    n = 101;
fxd = ParallelSum[(n/(2 \[Pi]))^-3 NIntegrate[
      8 Fn[x, y, z, r, 0.01], {r, -5, -1}, 
      Method -> {Automatic, "SymbolicProcessing" -> 0}], {x, \[Pi]/
     n, \[Pi], (2 \[Pi])/n}, {y, \[Pi]/n, \[Pi], (2 \[Pi])/
     n}, {z, \[Pi]/n, \[Pi], (2 \[Pi])/n}] // AbsoluteTiming;
SumInt = Table[
    fxd[[2]] + 
     ParallelSum[(n/(2 \[Pi]))^-3 NIntegrate[
        8 Fn[x, y, z, r, 0.01], {r, -1, R}, 
        Method -> {Automatic, "SymbolicProcessing" -> 0}], {x, \[Pi]/
       n, \[Pi], (2 \[Pi])/n}, {y, \[Pi]/n, \[Pi], (2 \[Pi])/
       n}, {z, \[Pi]/n, \[Pi], (2 \[Pi])/n}], {R, -1, 1.5, 0.05}] // 
Res = ParallelTable[{-1 + (j - 1) 0.05, SumInt[[2]][[j]]}, {j, 1, 
    50 + 1}];
ListLinePlot[{Res[[All, {1, 2}]]}, 
 PlotStyle -> {Red, Thickness[0.01]}, ImageSize -> 500]

enter image description here

  • $\begingroup$ Yours code dosen't give me a Plot? $\endgroup$ Aug 7, 2020 at 16:33
  • $\begingroup$ @MariuszIwaniuk, I added the part of the code to plot the results, check the updated version. $\endgroup$ Aug 7, 2020 at 17:50
  • $\begingroup$ Thanks for updated version :) $\endgroup$ Aug 7, 2020 at 18:01
  • $\begingroup$ Just to defend my honor, my plot has an offset since I defined the value at $-1$ to be zero, since NIntegrate was giving me inconsistent results for the integral from $-5$ to $-1$. Anyway, nice job :) $\endgroup$
    – Hausdorff
    Aug 7, 2020 at 23:22
  • $\begingroup$ @Hausdorff, you did a grate job, thanks to your detailed answer I was able to figure out this method. $\endgroup$ Aug 8, 2020 at 13:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.