Q:

NIntegrate for extremely oscillatory function

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,
Ella

A:

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:
NIntegrate[i0,
  {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:
NIntegrate[i0,
  {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}  *)

