I have a question about how to use NIntegrate. I need to integrate a slightly different function thousands of times, one very much like the one I have posted below. The functions are highly oscillatory so I usually use the Cuba Vegas package, but I've been having so much difficulty getting Vegas to work consistently on my institute's research cluster that I'm thinking of going back to functions inherent to Mathematica. Another similar question on here had the answer to use Method -> "LevinRule", but in this case it takes hours and doesn't seem to get anywhere. Otherwise I can maybe just up the PrecisionGoal, but it will start taking a very long time... Vegas evaluates this very fast, to a good precision, is there anything similar in Mathematica?

Could anyone point me towards a quick and accurate integration strategy for this type of function please?

a = Sqrt[29.7/53.77];

    NIntegrate[(((0.00025203302433875766 Exp[-(((Tan[x1]^2 + Tan[y1]^2)/
                5.8564 + (Tan[z1]^2)/1.9321)^(1/2))/
          a]) (Exp[-((((Tan[x1])^2 + (Tan[y1] - 0.4)^2)/
               5.8564 + ((Tan[z1])^2)/1.9321)^(1/2))/
         a]) (Exp[-((((Tan[x2] - 0.1)^2 + (Tan[y2] - 0.1)^2)/
               5.8564 + ((Tan[z2] - 0.1)^2)/1.9321)^(1/2))/
         a]) (Exp[-((((Tan[x2] - 0.1)^2 + (Tan[y2] - 0.1 - 0.4)^2)/
               5.8564 + ((Tan[z2] - 0.1)^2)/1.9321)^(1/2))/a]))/
    Sqrt[(Tan[x1] - Tan[x2] + 0.1)^2 + (Tan[y1] - Tan[y2] + 
         0.1)^2 + (Tan[z1] - Tan[z2] + 0.1)^2])*
  Sec[x1]^2 Sec[y1]^2 Sec[z1]^2 Sec[x2]^2 Sec[y2]^2 Sec[
    z2]^2, {x1, -Pi/2, Pi/2}, {y1, -Pi/2, Pi/2}, {z1, -Pi/2, 
  Pi/2}, {x2, -Pi/2, Pi/2}, {y2, -Pi/2, Pi/2}, {z2, -Pi/2, Pi/2}]

Here is an image of a couple of these integrations for similar functions (the x axis is just varying a constant so that we can see a couple of different results of the integration): Red - Vegas, Blue - MonteCarlo, Green - AdaptiveMonteCarlo, Black - MonteCarloRule


How do I know which of these is correct? AdaptiveMonteCarlo and MonteCarloRule both seem to be giving similar results which are lower than the similar results given by Vegas and MonteCarlo.

Thank you in advance.

Kind regards,


  • 2
    $\begingroup$ Can you show us the link for Cuba Vegas package? For this specific example, What's the result given by this package? How long does it take? $\endgroup$
    – xzczd
    Commented Jun 17, 2018 at 14:14
  • 3
    $\begingroup$ What does "slightly different" mean? $\endgroup$ Commented Jun 17, 2018 at 14:23
  • 1
    $\begingroup$ @xzczd Probably this: feynarts.de/cuba $\endgroup$
    – Michael E2
    Commented Jun 17, 2018 at 14:32
  • 3
    $\begingroup$ Cuba Vegas uses stochatic integration. So try out Method -> "MonteCarlo and friends. What is the precision you need for these integrals? $\endgroup$ Commented Jun 17, 2018 at 14:35
  • $\begingroup$ xzczd, yes the link posted by MichaelE2. @HenrikSchumacher, thank you for putting me on the right track - I'm trying the MonteCarlo methods now and they are matching the results that I had obtained with Vegas. However, how should I choose between AdaptiveMonteCarlo (gives result 0.0106408 for the function in the post), MonteCarlo (0.0110414) or MonteCarloRule (0.0109079)? I have added a picture comparing the different methods in the question - how do I know which method is giving the most real result? Thank you though, everybody. $\endgroup$
    – E Crane
    Commented Jun 17, 2018 at 16:47

1 Answer 1


The integrand is not oscillatory, but it is very high-dimensional. It does decay rather quickly at the boundary. (After an arctangent substitution, all the trig. functions go away and we have a exponentially decaying integral over ${\bf R}^6$.)

One can raise "MaxErrorIncreases" as suggested in the warning message or use Method -> "LocalAdaptive". By raising "MaxErrorIncreases", the integral converges: the error intervals are nested and decrease in length, but somewhat slowly. The "LocalAdaptive" method gives a result without warning messages but one that is somewhat different than the "GlobalAdaptive" method.

For a call of the following form, where i0 is the OP's integrand and "MaxErrorIncreases" -> 10000 may be set to whatever number, the results are summarized in the table below the call:

  {x1, -Pi/2, Pi/2}, {y1, -Pi/2, Pi/2}, {z1, -Pi/2, Pi/2},
  {x2, -Pi/2, Pi/2}, {y2, -Pi/2, Pi/2}, {z2, -Pi/2, Pi/2},
  Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> Incr}] // AbsoluteTiming

Incr    Integral     Error          Interval                    Timing
  2000  0.000807653  5.03281*10^-5  {0.000757325, 0.000857981}     9.31232
 10000  0.000795684  1.51190*10^-5  {0.000780565, 0.000810803}    49.7873
 50000  0.000792272  3.10038*10^-6  {0.000789172, 0.000795372}   259.574
200000  0.000791696  7.78874*10^-7  {0.000790917, 0.000792475}  1043.42

The "LocalAdaptive" result does not lie in all the above intervals:

  {x1, -Pi/2, Pi/2}, {y1, -Pi/2, Pi/2}, {z1, -Pi/2, Pi/2},
  {x2, -Pi/2, Pi/2}, {y2, -Pi/2, Pi/2}, {z2, -Pi/2, Pi/2},
  Method -> "LocalAdaptive"] // AbsoluteTiming
(*  {84.5794, 0.000806581}  *)

If the arc tangent substitution is made, I get a slight better result after a longer computation:

(*  {258.019903, 0.0007963143048912211}  *)
  • $\begingroup$ Thank you very much @MichaelE2. Why is Global or Local adaptive methods a better choice than MonteCarlo methods here? Also you seem to know which gives the best result by looking at the error, is there a simple way to display the error so that I can also compare strategies - in another answer you posted something to do with IntegrationMonitor? $\endgroup$
    – E Crane
    Commented Jun 17, 2018 at 18:15
  • 1
    $\begingroup$ @ECrane The error estimates are reported in the NIntegrate::eincr warning message, but you can also get it from IntegrationMonitor. I get messages in all four "GlobalAdaptive" cases, so reading the error off the messages was easiest. I don't have a good understanding of the strengths and weaknesses of Monte Carlo, except they are recommended for high-dimensional integrals because everything else often fails. Here we have a function that decays at the boundary, so I felt the standard multidimensional rule had a chance.... $\endgroup$
    – Michael E2
    Commented Jun 17, 2018 at 18:59
  • 1
    $\begingroup$ There is a singularity where the Sqrt in the denominator vanishes. You can point out the singular manifold to NIntegrate but the computed value of the integral seems stable. So I think that NIntegrate is handling the singularity satisfactorily. Admittedly I'm just relying on some heuristic evidence for this conclusion. $\endgroup$
    – Michael E2
    Commented Jun 17, 2018 at 19:02
  • $\begingroup$ thank you. That made it very clear. $\endgroup$
    – E Crane
    Commented Jun 18, 2018 at 19:45

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.