Earlier this year I reached an impasse with Maple's numerical integration package over the need for a high-resolution (say 50 or more digits) evaluation of multiple integrals for a Bayesian statistics application. An integral such as sqrt(1.+4.852800000000*chi^1411+5.887416960000*chi^200+4.852800000000*chi^100+5.887416960000*chi^2822-21.48043392000*chi^100*chi^1411) over a rectangular range like chi=0..0.9999,chi=0..0.9999999 would return unevaluated. Infolevel[`evalf/int`] explained that the methods available in Maple would not return the required precision. But the problem persisted even when I swallowed my pride and settled for 12 digits of precision. Stuck as I am, I haven't even tried to present the problem to Maple in higher dimensions.
To be clear: those were integer exponents under the radical, and the region of integration is rectangular but not square, between zero and almost one in each dimension. I'm also interested in the volume beneath this radical the other implied partitions of the unit square, i.e., (0.9999 to 1) x (0.999999 to 1) etc. I'm not merely summing these areas to together; they are of interest in their own right.
In general, the regions of integration are not related to the size of the exponents under the radical. I'm not religious about what methods I use in the various regions; mixing and matching is fine so long as they don't take forever and can provide some accuracy.
In an exchange this past March with Axel in this forum and in his blog, Fred Chapman mentioned the work that he and colleagues have done to develop the Geddes-Newton method for efficient multiple integration. I would love to try it. I see the papers, the thesis, and I have to think that there is a toolbox somewhere with a nice implementation of it, but I can't seem to find it. Can anyone point me in the right direction?
Thanks all (and congratulations to Fred and his colleagues on their contribution!),