Integration strategies for oscillatory multidimensional function

Question

I am seeking to integrate a highly oscillatory, multidimensional function. I am currently using NIntegrate's QuasiMonteCarlo approach. However, this is time-consuming and, given my current resources, not very accurate. How can I obtain more reliable estimates of the beasty integral given below? As the function itself will be integrated at a later stage, I am also interested in speeding up the function evaluation.

The integral to be solved:

fun[r_?NumericQ, d_?NumericQ, c_, opts:OptionsPattern[]]:=
   NIntegrate[
      Cos[(c (d^2+r^2-2 d r Cos[ta] - 3 ((-r+d Cos[ta]) Cos[tb]+d Cos[a] Sin[ta] Sin[tb])^2)) / Abs[d^2+r^2-2 d r Cos[ta]]^(5/2)] 
      Cos[(c (d^2+r^2+2 d r Cos[ta] - 3 ((r+d Cos[ta]) Cos[tb]+d Cos[a] Sin[ta] Sin[tb])^2)) / Abs[d^2+r^2+2 d r Cos[ta]]^(5/2)]
      * Sin[ta]*Sin[tb]/(2*Pi),
   {ta,0,Pi},
   {tb,0,Pi/2},
   {a,0,Pi},
   Evaluate@FilterRules[{opts},Options[NIntegrate]]
]

Typically, $d$ = 2, $r$ is in the range from 0 to Infinity (with the small values and $r$=$d$ posing problems), and $c$ is in the range from 100 to 3000. A typical function call is:

AbsoluteTiming[fun[3, 2, 400, Method -> "QuasiMonteCarlo", PrecisionGoal -> 6, 
MaxPoints -> 40000000]]
(* -> {102.3215798,-0.00442278} *)

This issues a NIntegrate::maxp warning and indicates an error estimate of 0.00011. Using the default strategy I obtain:

AbsoluteTiming[
 fun[3, 2, 400, MaxRecursion -> 20, Method -> {GlobalAdaptive, MaxErrorIncreases -> 10000}]]
(* -> {9.3912165,-0.00439357} *)

and a NIntegrate::eincr warning. Estimated error: 0.0369.

How to proceed from here? Thank you for your help.

Answer 1

Since no comment or answer has given as of the time I saw this post, and I don't have enough reputations to leave a comment, let me give a quick answer I used to solve the same kind of problems. I'm sure there must be a better way to do it with Mathematica, so this is just a beginning.

To evaluate an N-dimensional integral with a highly oscillatory integrand (which may peak at a particular small region in the N-dimensional space), I think all built-in methods and strategies of NIntegrate are not applicable. They are either super slow, or give a very inaccurate result which is easy to be proven wrong. Particularly, although to my knowledge Monte Carlo seems to be the best candidate for integrating a general multi-dimensional function (i.e., not knowing the symmetries in the integration region, so can't further simplify it and/or attack it with other strategies designed for special purposes), only when the integrand is smooth enough can Monte Carlo give a relatively good estimate and fast convergence. As a result, I would suggest using so-called "importance-sampling" Monte Carlo, which is not provided in built-in functions yet. (That's odd!)

Fortunately, Thomas Hahn (the author of FormCalc, FeynArts, etc.) has implemented this kind of algorithm in Mathematica (and C/C++ and Fortran as well). It's called the Cuba library. It provides both importance sampling and stratified sampling algorithms, and also the mixture of both. It is still being developed and under maintenance, so at least it works with Mathematica 8. Particularly, I would suggest using its Vegas routine, which is exactly designed for evaluating multi-dimensional integral, to solve your problem. Since it has a good user manual, I would simply refer you to that without further explanation.

However, if you don't mind doing the integral with C/C++, I would say that the Vegas routine provided in the GNU Scientific Library (GSL) converges way faster than the Cuba library does. Not sure why, maybe it has to do with the MathLink interface. Anyway, this is beyond the scope of this site.