hcubature http://ab-initio.mit.edu/wiki/index.php/Cubature_(Multi-dimensional_integration) seems to be a useful package for computing multidimensional integrals numerically. However, I have come upon several examples where hcubature fails. Consider the following contrived example:

$$\int_0^{10}dx\int_0^{10}dy\int_0^{10}dz\sin(x + y)\cos(y - z)\exp(-x*y^3)$$

When one computes this integral using hcubature, the (incorrect) result is 0.879338. However, when one splits the integral into two,

$$\int_0^{10}dx\int_0^{5}dy\int_0^{10}dz\sin(x + y)\cos(y - z)\exp(-x*y^3)\\+\int_0^{10}dx\int_5^{10}dy\int_0^{10}dz\sin(x + y)\cos(y - z)\exp(-x*y^3)$$

and then evaluates each term using a call to hcubature, the result is 0.892735, and I presume this to be the correct result.

One can see the problem here by plotting the points at which hcubature evaluates the integrand. For the single hcubature calculation on the entire 10x10x10 domain, the evaluation points projected on the x-y plane are plotted here. enter image description here

hcubature fails since the adaptive scheme does not place any evaluation points in the region $5<y<10,0<x<0.2$.

While the correct answer of 0.892735 is also achieved by spicy.integrate.tplquad, and also by Mathematica with the following call

 Sin[x + y]*Cos[y - z]*Exp[-x*y^3], {x, 0, 10}, {y, 0, 10}, {z, 0, 
  10}, Method -> "GaussKronrodRule"]

the problem is not unique to hcubature. Using Monte Carlo integration in Mathematica gives us a result of 0.878907 with the following call

 Sin[x + y]*Cos[y - z]*Exp[-x*y^3], {x, 0, 10}, {y, 0, 10}, {z, 0, 
  10}, Method -> "MonteCarloRule"]

and presumably the MonteCarloRule is failing for the same reason as hcubature.

I presume that the failure to miss specific regions of the domain where the integrand is nontrivial is a general challenge for adaptive methods, even for analytic functions like the one I have proposed. Then how can one be confident when using an adaptive quadrature method that the result is accurate, especially when one is doing thousands of integrals and cannot visually inspect the sampling of each one.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.