Cubature formulas for the unit disk

This page is about cubature formulae for numerical integration over the unit disk. In the first part we try to derive coefficients by direct approach, second is devoted to Gauss rules for 2D unit disk based on product of two 1D quadratures.
We offer C/C++ open source library which allows user to estimate integral over disk $x^2+y^2\le R^2$ with center at $(Xc,Yc)$ where $R, Xc, Yc$ can be arbitrary. Algorithm is based on Gauss product rules and includes high-precision coefficients.

Direct derivation

We would like to build approximation of the integral

$I(f) = \displaystyle{\iint\limits_{x^2+y^2\le 1}\,f(x,y)\,dxdy}$

by the weighted sum of integrand values

$C(f) = \sum_{i=1}^N\omega_i\,f(x_i,y_i)$

where $\{\omega_i,x_i,y_i\}$ are such that $C(f)$ is exact on polynomials $x^iy^j,\,\,0\le i+j\le n$ . Formulas of such type are called cubatures of algebraic degree $n$ .

After switching to polar coordinates we obtain

$I(f) = \displaystyle{\iint\limits_{x^2+y^2\le 1}\,f(x,y)\,dxdy}=\int_0^1\int_0^{2\pi}\,\rho\,f(\rho\cos\varphi,\rho\sin\varphi)\,d\varphi d\rho$

and

$C(f) = \sum_{i=1}^N\omega_i\,f(\rho_i\cos\varphi_i,\rho_i\sin\varphi_i)$

Let’s represent $f(x,y)$ in the form of multivariate Maclaurin series as

$\begin{aligned} f(x,y) & \approx f(0,0) +x\, f_x(0,0) +y\, f_y(0,0) \\ & {}\quad + \frac{1}{2!}\left[ x^2\,f_{xx}(0,0) + 2xy\,f_{xy}(0,0) +y^2\, f_{yy}(0,0) \right]+\cdots, \end{aligned}$

This gives us approximation error $Z(f) = I(f)-C(f)$ conveniently expressed through function derivatives.

To unsure exactness on polynomials up to total degree of $n$ we have to find $\{\omega_i,\rho_i,\varphi_i\}$ such that multipliers of all partial derivatives $\displaystyle\frac{\partial^{i+j}f}{\partial^ix\partial^jy },\,\,0\le i+j\le n$ in $Z(f)$ are equal to zero. This leads to the system of nonlinear equations, which can be solved analytically for small $n$ and numerically for the higher degrees.

Additionally we set $\{\varphi_i = \frac{2\pi}{N}(i-1)\}$ and $\{r_i\}$ to be symmetric, since such choice makes cubature exact on $x^iy^j$ naturally, where at least one of $i$ or $j$ is odd.

Below we show several cubature formulas derived by the explained procedure.

$\renewcommand*\arraystretch{2.5} \begin{tabular}{|r|c|c|} \hline $N$& $N$ point cubature $\{\omega_i,\rho_i,\varphi_i = \frac{2\pi}{N}(i-1)\}$ & Approx. Error \\ \hline 4&$\displaystyle{\rho_i=\frac{\sqrt{2}}{2},\,\,\omega_i =\frac{\pi}{4},\,\,i=1\hdots N}$ & $O(f^{(4)})$\\ \hline 8&$\displaystyle{\rho_{2k-1}=\sqrt{\frac{3+\sqrt{3}}{3}},\,\,\omega_{2k-1}=\frac{\pi(2-\sqrt{3})}{16}}$&\\ &$\displaystyle{\rho_{2k}=\sqrt{\frac{3-\sqrt{3}}{3}},\,\,\omega_{2k}=\frac{\pi(2+\sqrt{3})}{16}}$&$O(f^{(6)})$\\ &$k=1\hdots N/2$&\\ \hline 16&$\displaystyle{\rho_{2k-1}=\sqrt{\frac{3+\sqrt{3}}{6}},\,\,k=1\hdots N/2}$&\\ &$\displaystyle{\rho_{2k}=\sqrt{\frac{3-\sqrt{3}}{6}},\,\,k=1\hdots N/2}$&$O(f^{(8)})$\\ &$\displaystyle{\omega_i=\frac{\pi}{16},\,\,i=1\hdots N}$&\\ \hline \end{tabular}$

Please leave comment on this page if you find these formulas useful.

Product rule

Derivation of optimal cubatures with minimal function evaluations required to achieve given algebraic degree is not an easy task, they are known only for some $N=3,4,6,7,10,12,18$ (see [1],[2]).

However if we want to build cubatures for any desired approximation order we can use product of two 1D rules to estimate $I(f)$ . Although this approach is not optimal in the sense of number of integrand evaluations it gives flexible and easy way of cubature derivation (see [1],[2],[6] for details) which we repeat here.

At first we use slightly different transformation to polar coordinates:

$I(x^i\,y^j) = \displaystyle{\iint\limits_{x^2+y^2\le 1}\,x^i\,y^j\,dxdy}=\int_{-1}^{1}\,|\rho|\rho^{i+j}d\rho\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\,{\cos^i\varphi}\,{\sin^j\varphi}\,d\varphi$

After substitution $t=\sin\varphi$ in the second integral we get

$I(x^i\,y^j)=\displaystyle{\int_{-1}^{1}\,|\rho|\rho^{i+j}d\rho\int_{-1}^{1}\,(1-t^2)^{-1/2}\,(1-t^2)^{i/2}\,t^j\,dt}$

These integrals can be computed separately from each other by using 1D Gauss quadratures.

We use Gauss-Chebyshev formula of the first kind for approximation

$\displaystyle{\int_{-1}^{1}(1-t^2)^{-1/2}\,f(t)\,dt \approx \sum_{i=1}^M B_i\,f(t_i)}$

where nodes and weights are known in the closed form

$\displaystyle{t_i=\cos\left(\frac{2i-1}{2M}\right)\pi\qquad B_i=\frac{\pi}{M},\qquad i=1,\hdots,M}$

To compute second integral

$\displaystyle{ \int_{-1}^{1}\,|\rho|\,f(\rho)\,d\rho \approx \sum_{j=1}^M A_j\,f(\rho_j) }$

we need to build Gauss quadrature with weight function $|\rho|$ . I have no access to [3] where it is considered in details, so I derived points and coefficients $\{A_j,\rho_j\}$ using standard procedure for constructing quadrature formulas [5]:

Find coefficients of two-term recursion relation for monic polynomials $P_n(x)$ orthogonal on $[-1,1]$ with weight $|x|$ :
$\displaystyle{P_{n+1}(x)=(x-a_n)P_n(x)-b_n\,P_{n-1}(x)}$
Find eigenvalues and first component of eigenvectors of $M\times M$ principal leading minor of Jacobi matrix:
$\displaystyle{ \left[\begin{matrix} a_0 &\sqrt{b_1}&&&&0\\ \sqrt{b_1}&a_1 &\sqrt{b_2}&&&\\ &\sqrt{b_2}&a_2 &\sqrt{b_3}&&\\ 0 & &\ddots &\ddots &\ddots& \end{matrix}\right] }$
It is a symmetric tridiagonal matrix. Implicit QL algorithm can be applied to solve this task efficiently (sterf, stein routines in LAPACK). Then points $\{\rho_j\}$ are calculated eigenvalues. Weights $A_j=\mu_0(\phi_1^j)^2$ where $\phi_1^j$ is the first component of corresponding eigenvector).

In our case $a_n=0,\,\,\mu_0=1$ , and $b_n$ can be calculated by analytical formula for any $n$ (send me e-mail if you are interested).

Combining both rules $\{A_j,\rho_j\}$ and $\{B_i,t_i\}$ we obtain Gauss-type cubature (or product rule) for the unit disk which has $N=M^2$ points and which is exact on polynomials of a total degree at most $2M-1$ . Distribution of the points for some $M$ can be observed on these plots (click on image to zoom):

Download

C/C++ library implementing Gauss product rule for the integration over the disk (2D sphere) is available for download:

Gauss 2D Sphere Cubature.zip (400KB)

It uses high-precision points and coefficients of 1D rules $\{A_j,\rho_j\}$ and $\{B_i,t_i\}$ to approximate integral over 2D sphere $x^2+y^2\le R^2$ with center at $(Xc,Yc)$ .

Tools

I have used Maple for all symbolic calculations and QuickLaTeX for rendering mathematical equations on the web page. Unfortunately Maple was too slow and unstable for derivation of 1D $|\rho|$ -quadrature with $M>20$ . I ended up writing my own C++ program specifically for that using ALGLIB and multi-precision arithmetic library MPFR. Let me know if your interested in this program.

References

The most complete and recent compendium of cubature formulas is Encyclopaedia of Cubature Formulas maintained by Ronald Cools . He wrote several surveys on cubature formulas, some of them are available online, check his publication page.

In particular his survey on cubatures for unit disk is available here A survey of known and new cubature formulas for the unit disk. Although it doesn’t contain numeric values for weights and points (exept for $N=57,69,88$ ) it has huge number of useful links to other topic-specific papers.

[1] A.H. Stroud, Approximate calculation of multiple integrals, 1971.
[2] R. Cools, K.J. Kim, A survey of known and new cubature formulas for the unit disk.
[3] A.H. Stroud, D. Secrest, Gaussian Quadrature Formulas, 1966.
[4] G.H. Golub, J.H. Welsch, Calculation of Gauss quadrature rules, Math. Соmр. 23 1969.
[5] W. Gautschi, Orthogonal Polynomials: Computation and Approximation, 2004.
[6] I.P. Mysovskikh. Interpolatory cubature formulas. Moskva: “Nauka”. 1981. (Russian).

(8 votes, average: 4.50)

13 Comments

Nagaraju Krishnappa
Posted October 24, 2011 at 10:51 pm | #
Hi
I downloaded Gauss_2D_Sphere_Cubature.zip and try to run example.c but got an error message that gauss_product_2D_sphere is not found.
Could you please let me know if gauss_product_2D_sphere is missing in the distribution and if it is send set to me.
Thanks
nagaraj.
- Pavel Holoborodko
  Posted October 25, 2011 at 9:51 am | #
  You have to compile & link gauss_2d_sphere.c along with example.c.
Brian
Posted December 19, 2012 at 1:57 am | #
Thanks for making this available. When running the example, I notice that the error in the integral starts to increase for n > 50 (10^-13 as opposed to 10^-15). I don’t need that level of precision, but was curious if this is something particular to the function used in the example or a generic effect. Also, one rule in particular stands out: “n = 384: error = 2.44418289875682”. Maybe one of the coefficients for this rule is incorrect?
- Pavel Holoborodko
  Posted December 19, 2012 at 11:14 am | #
  Brian, thank you for your feedback! My code contains cubature of following orders:
  1-50, 64, 128, 192, 256, 320, 448, 512, 576, 640, 704, 768, 832, 896, 960, 1024, 2048.
  No 384 is supported – sorry for not making clear this on the page.
  - Brian
    Posted December 19, 2012 at 5:06 pm | #
    Ahh. Sorry, I assumed since it was being run in the example that a rule existed. The fact that the error is larger than the true value should have been enough of a hint…
Bernard SIRET
Posted February 23, 2013 at 5:45 am | #
Hi
I am trying to develop gauss like formulas for cubature over the unit disc with the constraint that all the points need to be located on three straight lines, all three lines intersecting on the boundary; in one point (say x=0, y=-1). This leaves entire areas not covered but what can be done?
- Pavel Holoborodko
  Posted February 23, 2013 at 9:35 pm | #
  Hi Bernard,
  Interesting case. You can build 2D interpolating polynomial over the points you have and then derive cubature formula by integrating this polynomial. Take a closer look to “Direct derivation” section on the page – it briefly describes this idea (using power series instead of polynomial).
  Let me know about your further progress on the task – it is very interesting for me.
  By the way, are there any restrictions on number of points?
  Location of lines is fixed or can be chosen?
  Location of points on the lines is fixed (uniform sampling) or can be chosen?
  Let me know if you need my further help.
Bernard SIRET
Posted February 25, 2013 at 12:46 am | #
This is a practical industrial problem
We want to get a flow by “sampling” velocities, in a pipe, in cases the ISO standart cannot be used
So practically the number of points is less than 16
The position location of the lines can be free as long as all 3 lines intersect on the circle x2+y2=1
more than 3 lines could be considered, all n lines having a common intersect, on the unit circle
On each line the location of the points is free. All points must be inside the unit disk, strictly (x2+y2 < 1)
The interpolating polynomial is not practical for field people. Setting a probe at given, preset locations, &nd appplying weights is. Regards
John Pennucci
Posted August 21, 2013 at 5:42 am | #
In your example for an 8 point cubature the radius of 4 of the points is sqrt((3+sqrt3)/3). This is obviously greater than 1.0 and outside the unit circle.
- Pavel Holoborodko
  Posted September 10, 2013 at 10:39 pm | #
  It is perfectly normal for cubature to have nodes outside unit disk. Check out paper by R. Cools (second reference in the list).
  Constructing of cubature with all nodes inside sphere – is a real challenge and might not be possible at all (at least for algebraic accuracy)
Jason Williams
Posted March 12, 2017 at 2:31 pm | #
Pavel,
Many thanks for freely providing the information on this webpage. I would like to use this cubature in my meshfree code and would like to confirm the points and weights. Are the cubature points given by x = r[j]*q[i] and y = r[j]*t[i] as shown in the example C++ code? With that said are the weights given by A*B? I am also interested in the C++ program using ALGIB mentioned above.
Many thanks,
JW.
marco
Posted March 19, 2017 at 6:16 am | #
Hi,
thanks for the article! I have a question: what if I want to integrate over a regular multidimensional rectangular domain, let’s say
$X = [x_{1,lower} , x_{1,upper}] \times [x (x_{2,lower} , x_{2,upper}] \times \dots\times [(x_{n,lower} , x_{n,upper}]$
and I want to compute the following definite n-dmensional integral
$\int_{x_{1,lower}}^{x_{1,upper}}\int_{x_{2,lower}}^{x_{2,upper}} \dots \int_{x_{n,lower}}^{x_{n,upper}} f(x_1,x_2,\dots,x_n)dx_1 dx_2\dots dx_n$
what’s the best numerical approach known in literature?
thanks!
marco
Nick
Posted December 16, 2017 at 6:47 am | #
Hi Pavel,
Great resource! Question: Do you know where I could find Kronrod-type extensions for these formulae? Cools lists a couple of rules in his survey, but they all have fatal flaws, like negative weights, nodes outside the unit disk or on the boundary. Do Kronrod style extensions with positive weights and interior nodes exist?