# Numeric integration of multi-dimensional integral with known boundaries

I have a (2-dimensional) improper integral

$$I=\int_A \frac{W(x,y)}{F(x,y)}\,\mbox{d}x\mbox{d}y$$

where the domain of integration $A$ is smaller than $x=[-1,1]$, $y=[-1,1]$ but further restricted by $F(x,y)>0$ . Since $F$ and $W$ are smooth and $W \ne 0$ at the boundaries, the later relation implies that the integrand can be singular at the boundaries. The integrand is finite though. I, so far, compute this integral with nested numerical integration. This is successfull but slow. I search a more appropriate (faster) method to address the integral, maybe a Monte-Carlo method. But I need one which does not put points on the boundary of the non-cubic domain A and takes the limit of the improper integral correctly. Can a Integral transformation help for this general expression? Note that I can solve $F(x,y)$ for $y$ as a function of $x$ and even compute $I$ for a few special weight functions $W(x,y)$.

-
 Can you be a bit more specific about what methods you've used so far? What specific routines have you used in a nested way? Also, is $F(x,y) \ne 0$ inside $A$, i.e. are the roots of $F(x,y)$ only on the boundary? – Pedro Aug 23 '12 at 8:38 The GSL algorithm QAGS: gnu.org/software/gsl/manual/html_node/…. Thanks for the edits (didn't see the option to typeset equations)! – highsciguy Aug 23 '12 at 8:51

Disclaimer: I wrote my PhD thesis on adaptive quadrature, so this answer will be severely biased towards my own work.

GSL's QAGS is the old QUADPACK integrator, and it is not entirely robust, especially in the presence of singularities. This usually leads to users requesting far more digits of accuracy than they actually need, thus making the integration quite expensive.

If you're using GSL, you might want to try my own code, CQUAD, described in this paper. It is designed to cope with singularities, both at the interval edges and within the domain. Note that the error estimate is quite robust, so only ask for as many digits as you actually require.

With regards to Monte-Carlo integration, it depends on what kind of accuracy you're looking for. I'm also not quite sure of how well it will work near singularities.

-
I will definitely have a look at this because it will be most simple to implement it. In fact I experienced that the QAGS routine was not super-stable for this problem. – highsciguy Aug 23 '12 at 9:55
Is there a way to influence the occurance of 'GSL_EDIVERGE'? It seems to appear for some parameters. – highsciguy Aug 23 '12 at 12:31
@highsciguy: The algorithm returns GSL_EDIVERGE when it believes that the integral is not finite. If you could give me an example for which it fails, I could have a closer look at it. – Pedro Aug 23 '12 at 13:15
Its a little difficult to isolate a simple routine, as it is embeded in a generic code for n-dimensinal integrals. I will see ... But for fixed y, 1/sqrt(F(x,y)) should behave like 1/sqrt(x) as x approaches the zeros of F(x,y) since F(x,y) can then be written as a polynomial in x. But could be that the 1/sqrt(x) behaviour starts late. Could also be that the numeric precission of the integrand is not too good. – highsciguy Aug 23 '12 at 16:44
@highsciguy: Yes, this is a bad idea. Most quadrature rules assume that the integrand has some degree of smoothness, and if you set it to zero as of some arbitrary point, you're introducing a discontinuity. You will get much better results if you use the actual interval! – Pedro Aug 24 '12 at 8:34
The reason is relatively easy to understand. Take, for example, just $\int_{[0,1]^n} f(x)\; d^nx$ where $n$ is the dimension of the problem. Let's say, for simplicity, that you subdivide every dimension into $M$ sub-intervals, i.e., you get $M^n$ hypercube cells in total. Let's assume further that you use a Gauss formula with $k$ Gauss points, just as an example. Then you have $N=(kM)^n$ quadrature points in total, and because the $k$ Gauss points provide you with order $(2k-1)$ accuracy, $e = {\cal O}(h^5)={\cal O}(M^{-(2k-1)})$, your overall accuracy as a function of evaluation points will be $$e = {\cal O}(N^{-(2k-1)/n}).$$ On the other hand, Monte Carlo methods generally only provide you with error convergence as $$e = {\cal O}(N^{-1/2})$$ which is worse than for any Gauss formula with at least $k>n/4+1/2$ points per interval. The reason is relatively simple to understand: Gauss quadrature chooses the interpolation points in some sort of smart way, Monte Carlo isn't. You can't expect anything useful to come of the latter. (Of course there are situations where Gaussian quadrature is difficult: for example in your case where the integration domain is irregularly shaped; but in that case you're likely still better off doing adaptive integration or similar.)
Now, there are practical (stability) problems with integration with more than, say, 8 or 10 points per interval. So if you want $k\le 8$, then you can't go beyond $n=30$. On the other hand, in that case, even choosing a single interval per direction ($M=1$) yields $N=8^{30}$ integration points, far more than you could ever in a lifetime evaluate. In other words, as long as you can evaluate enough integration points, quadrature on subdivisions of your integration domain is always the more efficient approach. It's cases where you have have a high dimensional integral for which you can't evaluate the integration points on even a single subdivision any more that people use Monte Carlo methods despite their worse convergence order.