Integration strategies for oscillatory multidimensional function

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.

-

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.

-
Well, NIntegrate[] does support stratified sampling (via Method -> "AdaptiveMonteCarlo"), but you are right that the lack of support for importance sampling is an unfortunate omission. Maybe future versions will have it... –  Ｊ. Ｍ. Jun 8 '13 at 16:48
You're always right @0x4A4D :) –  Leo Fang Jun 8 '13 at 16:51
@Leo Thank you for your recommendations. I second your observation regarding the issues of the built-in methods with oscillatory multidimensional integrals. I have eventually settled for Vegas as provided by the GSL, which I call through LibraryLink. Unfortunately, I have not yet been able to compile the Cuba library on windows. Using MSYS and the mingw compiler (64 bit) the make process bails out on linking (undefined reference to ´mkstemp´). Apparently, this is a known problem with mingw. Anyway, since the GSL approach lives up to my expectations, I did not try to resolve this issue ... –  user7935 Jun 10 '13 at 8:18
@Leo By the way, using 5000000 samples (1M for warm up) Vegas from the GSL yields -0.00442672 with error 0.0000624632 in 8.8 s. –  user7935 Jun 10 '13 at 8:21
@user7935, glad to know that it works for you! I don't use Windows here so I can't test the compilation...For the GSL part I roughly estimated it before, and I found that if numbers you expect to get are below ~10^-4 (so very close to zero), it may not be able to give a good estimate. Fortunately this seems not to be your case. BTW, just curious, are there really two cosine functions in the numerator of your integrand? –  Leo Fang Jun 10 '13 at 13:38