I'm trying to numerically integrate a multi-dimensional expression. The integrand is complicated; for example this is the integrand for $N=4$: $$\begin{aligned}&x_1^6x_2^5x_3^3x_4^2(x_1-x_1x_2)(x_1-x_1x_2x_3)(x_1-x_1x_2x_3x_4)\times{}\\&{}\times(x_1x_2-x_1x_2x_3)(x_1x_2-x_1x_2x_3x_4)(x_1x_2x_3-x_1x_2x_3x_4)\end{aligned}$$

And I'm most interested in expressions that have even larger $N$. The integration limits in my trial case are $x_1$ from 0 to 10, and everything else from 0 to 1.

To numerically integrate this I'm using Mathematica while double checking with the VEGAS integrator in C++ (via the GNU scientific library). The problem is the two integrators match for lower dimensions but not for higher ones. The agreement is great for $N =2, 3$, but for $N=4, 5$ there's a considerable discrepancy of about a factor of 5 that gets worse with higher $N$. Agreement doesn't fail completely; the results are still of the same order of magnitude, but a factor of 5 is too large for comfort especially with both integrators reporting error bars much lower than this.

I am unfamiliar with the internals of both integrators and how numerical integration in higher dimensions works. Is there a (non-bug) reason for why the two integrators give the same results at low dimensions but not at higher dimensions? If the only possible explanation is "a bug in your code" that's fine, but if there's some other possibility I'd like to check that out too.

  • 3
    $\begingroup$ The first thing that comes to mind is that you have very high powers of $x$, so over $[0,1]$ (assuming that's the interval) it's a very lopsided integrand, making it more difficult for Monte Carlo methods. So try a substitution that gets rid of pure powers of variables, i.e. $x_1=u_1^{1/13}$, etc., then maybe factor out the very small constant factor you'll get because that can mess with default tolerances, (like this), and check if that does anything. $\endgroup$ – Kirill Jan 25 '18 at 3:23
  • $\begingroup$ What does your boundary look like? $\endgroup$ – user14717 Jan 25 '18 at 4:00
  • $\begingroup$ @user14717 $x_1$ from 0 to infinity, the rest are from 0 to 1. $\endgroup$ – Allure Jan 25 '18 at 6:33
  • $\begingroup$ @user14717 sorry, it should be $x_1$ from 0 to 10. The integral to infinity is because I've already excluded an exponential term that should be in the integrand as well. $\endgroup$ – Allure Jan 25 '18 at 6:41
  • 2
    $\begingroup$ Can't you just do this integral in closed form? It's just a long sum of products of powers. $\endgroup$ – Kirill Jan 25 '18 at 17:50

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.