# Numerical integration over 2D disk

I have a real-valued function $f$ on the unit disk $D$ that is fairly well behaved (real-analytic everywhere) and would like to find the integral $\int_D f(x,y)dxdy$ numerically. After much searching, the best methods I've been able to find have been the simple quadrature rules in Abramowitz and Stegun that sample $f$ at up to 21 points. What work has been done since? In particular, I'm interested in rules that allow sampling at more than 21 points. One reference indicates that finding optimal quadrature rules is a hard problem, but it seems to me that something better must have been published in the last 50 years.

A couple of references suggest integrating over various domains by triangulating them and using numerical integrals over the triangles. Is this the preferred method for a disk?

(I'm trying to improve code that has been implemented using quasi-monte carlo. It seems to me that we could do much better using the knowledge that $f$ is real-analytic and probably well approximated by polynomials.)

Update: I can't easily say exactly what the functions $f$ are as they're the messy result of a chain of computations. I can say that qualitatively it's like a gaussian with a central hump, fast decay, though not exact rotational symmetry. I do have a pretty good handle on how big the hump is and where it's centred. All variations on this might happen and some are easy to dispense with: eg. the hump might be situated well outside the disk so I know the integral is nearly zero. Or the hump may be very wide in which case the integrand is almost constant. Sometimes the hump is contained well within the disk in which case I can switch to more efficient quadrature over the (approximate) support of $f$ rather than the disk. But having said all that, I'd still like to see some general gaussian quadratures rules for the disk that would apply to integrating any function over the disk that is well approximated by a polynomial.

Update2: After much web searching I found some Fortran code to do what I want (and more) and a reference to a book by Arthur Stroud, Approximate Calculation of Multiple Integrals. It seems as that this work from 1971 is the state of the art.

• Might help to state what the function actually is... – Noldorin Jun 18 '10 at 21:50
• Am I missing something or why not just change variables from the disk to the unit square and use nested Gaussian quadrature? – Fredrik Johansson Jun 18 '10 at 22:10
• @Fredrik You mean use polar coordinates to map from the unit square to disk? That gives an inefficient way to sample the space as you end up with many more points near the centre of the disk. – Dan Piponi Jun 18 '10 at 22:46

See the Encyclopedia of Cubature Formulas. The site is password protected, but the maintainer will give a password to anyone who asks.

• That looks like exactly what I want. Password is a bit weird and I worry there'll be some IP issues. I'll see... – Dan Piponi Jun 19 '10 at 2:27
• The site doesn't seem to be password-protected anymore. – gmvh Mar 13 at 15:27

Here's another one. a little strange, since the quadrature rule involve the value of line integrals, but I guess you can write those using a separate quadrature rule...

Numerical integration over a disc. A new Gaussian quadrature formula

Borislav Bojanov and Guergana Petrova

" We construct a quadrature formula for integration on the unit disc which is based on line integrals over distinct chords in the disc and integrates exactly all polynomials in two variables of total degree . "

Cheers.

You might like this one, I do not know if the entire link will fit.

http://www.math.tamu.edu/~gpetrova/

Journal of Approximation Theory Volume 104, Issue 1, (May 2000)

Uniqueness of the Gaussian Quadrature for a Ball Pages 21-44

Borislav Bojanov and Guergana Petrova

Department of Mathematics, University of Sofia, Boulevard James Boucher 5, 1164, Sofia, Bulgariaf1

Department of Mathematics, University of South Carolina, Columbia, South Carolina, 29208, U.S.A., f2 Received 8 June 1999; accepted 22 October 1999. ; Available online 26 March 2002.

Abstract

We construct a formula for numerical integration of functions over the unit ball in Image d that uses n Radon projections of these functions and is exact for all algebraic polynomials in Image d of degree 2n−1. This is the highest algebraic degree of precision that could be achieved by an n term integration rule of this kind. We prove the uniqueness of this quadrature. In particular, we present a quadrature formula for a disk that is based on line integrals over n chords and integrates exactly all bivariate polynomials of degree 2n−1.

Author Keywords: Gauss quadrature formula; orthogonal polynomials; highest degree of precision

Different article by same people:

http://www.math.tamu.edu/~gpetrova/CAM7238.pdf

• I should have said in my original question "A web search quickly turned up a paper by Bojanov and Petrova but it involves line integrals which is a bit weird and not useful for me" :-) It's sort of in the right spirit though in the sense that it's looking at schemes to squeeze out integrals for the maximum degree of polynomial. – Dan Piponi Jun 18 '10 at 23:36
• Fair enough. She has a million papers on her web page, publications section, many downloadable (such as Yossi's answer). I would guess that the reference sections of some of these papers would be worthwhile. I cannot say I have ever programmed numerical integration. Have you looked at the standard cambridge.org/us/catalogue/catalogue.asp?isbn=9780521880688 – Will Jagy Jun 18 '10 at 23:46

Here is link to my site, where I publish some of the cubatures of such kind:

Cubature formulas for the unit disk

Page also gives C source code for the numerical integration by the product of two 1D Gauss-type quadratures. Images of points distribution in that case are available too.

Let me know if you find it useful.

• That is useful. But in a way it fails to do what I'd really like. Given a limited CPU time budget for sampling the integrand, I'd like to use samples as efficiently as possible. Your schemes, being product schemes, put much more effort into the centre of the disk than round the edge, and suffer from some redundancy. The best schemes should distribute samples more evenly and so extract more independent information per sample. – Dan Piponi Sep 6 '10 at 19:20
• Besides product rule description & source code page lists optimal cubatures for $N=4,8$ you are looking for (see table in <i>Direct derivation</i>). – Pavel Holoborodko Sep 7 '10 at 0:48

An insufficiently well-known (so, perhaps slightly beyond the state of the art) integration algorithm can be found in the paper of O. Jenkinson and M. Pollicott entitled "A dynamical approach to accelerating numerical integration with equidistributed points". They claim a qualitative improvement over the state of the art as of three years ago.

An interesting approach can be found in the paper Extensions of Gauss Quadrature Via Linear Programming by Ernest Ryu and Stephen Boyd.

Their method is to show that ordinary Gauss quadrature on the real line can be interpreted as a solution to a linear programming (LP) problem and then use the same LP method on domains other than the line. When used in a more general setting it's no longer solving precisely the same problem that Gauss quadrature is solving but it's still a good heuristic for generating efficient quadrature methods.

One nice thing about their approach is that it generalizes to irregular domains too.

Forgive me for stating the obvious, but your best bet will be to use matlab/octave/scilab/.... I'd assume that any modern general purpose code will be adaptive (subdivide the area into squares recursively, estimate the integral on each square using a simple quadrature rule and spend most effort on the areas that are converging most slowly).

I only know enough to be dangerous. I would not attempt to write my own code for this purpose, since the code produced by a real numerical analyst will be vastly superior.

• I'm trying to write code that runs fast and performs many of these integrals. Matlab etc. won't do at all. And I know from Abramowitz and Stegun that people have worked on schemes specifically for integrating over disks. – Dan Piponi Jun 18 '10 at 22:49
• "fast" as in "in the inner loop"? I would agree with Chris -- matlab/octave/scilab/mathematica/etc have reasonably state of the art algorithms (which are notoriously hard to implement, so I am glad someone is paid to do this). – Igor Rivin Dec 16 '10 at 3:09