12

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)$.

Pedro
  • 9,573
  • 1
  • 36
  • 45
highsciguy
  • 1,119
  • 9
  • 16
  • 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 08:38
  • The GSL algorithm QAGS: http://www.gnu.org/software/gsl/manual/html_node/QAGS-adaptive-integration-with-singularities.html. Thanks for the edits (didn't see the option to typeset equations)! – highsciguy Aug 23 '12 at 08:51

3 Answers3

8

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.

Pedro
  • 9,573
  • 1
  • 36
  • 45
  • 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 09: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
  • Ah, yes. Right now I am integrating over the whole domain [-1..1],[-1..] and set the integrand zero outside. Is this an issue? – highsciguy Aug 23 '12 at 16:52
  • 1
    @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 08:34
  • I am approaching a sulution where the singularities of the integrand are computed numerically and the integral is split into multiple integrals between the singularity. It would be good to know precise the determination of the singularities needs to be. I.e. can the finite precision in the determination of the singularities be dangerous because the integration routine puts grid points too close to the singularities so that the numerical error in their determination matters? – highsciguy Sep 08 '12 at 22:16
  • @highsciguy: That depends on how the code deals with non-numerical, e.g. NaN and Inf, function values. My own code removes such nodes from the interpolant/quadrature rule. Matlab's routines simply fail and a number of other routines replace these values with zeros. In most cases, it should suffice to be close enough such that the node on the edge of the domain with the singularity evaluates to NaN or Inf and let the quadrature routine, whichever you're using, do its thing. – Pedro Sep 09 '12 at 08:10
  • I am doing some debuging now on my code for the nested integral. So far the CQUAD routine behaves worse (the results are less accurate than those obtained by QAGS with identical arguments). I would like to understand this. Hence the following question: Is possible to read and output the function evaluation abscissas and results from the different GSL workspace structures. I would like to see where the different algorithms evaluate the integrand. I tried to find it using a debugger and by looking into the code but could not access the arrays. – highsciguy Sep 13 '12 at 09:16
  • @highsciguy: You could wrap your integrand in another function which stores and/or outputs the nodes at which it is called. This is probably the easiest way to access this data. As for the internals, there are some debugging printfs in the source code that you could uncomment. Could you send me an example for which CQUAD does worse than QAGS? I'd be very interested in seeing how and why it does worse. – Pedro Sep 13 '12 at 09:29
  • Ah, yes I forgot to mention that, for the previous question, I was assuming that the integrand is set to zero in any case if the result would be NaN or Inf otherwise. But this means that the integrator could see a step in the interval if it looks too close by the singularities. But your answer means probably that it would be better to return NaN for your algorithm. – highsciguy Sep 13 '12 at 09:29
  • @highsciguy: Yes, the algorithm can deal with NaNs and Inf, so just let your code do its thing. Are you also integrating only in the interval of interest, e.g. you are no longer setting the values outside the actual interval to zero? – Pedro Sep 13 '12 at 09:31
  • I will try to extract such a case from the multiple integral. But right now there is still a possibility that I am making a mistake ... I will follow your proposal, just thought that it might be simpler to read it from the workspaces. – highsciguy Sep 13 '12 at 09:34
  • I integrate only inside the interval. – highsciguy Sep 13 '12 at 09:35
5

Monte Carlo methods can in general not compete with adaptive quadrature unless you have a high dimensional integral where you cannot afford the combinatorial explosion of quadrature points with the dimension.

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.

Geoff Oxberry
  • 30,394
  • 9
  • 64
  • 127
Wolfgang Bangerth
  • 55,373
  • 59
  • 119
1

Try a nested Double-Exponential Quadrature (see the implementations of Ooura). This technique use a variable transformation that makes the transformed integrand to behave very smoothly at the boundaries and is very efficient for handling singularities at the boundary. There is also a very good list of references on the DE quadrature on his website.

GertVdE
  • 6,179
  • 1
  • 21
  • 36